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SUMMARY 


In  this  annual  report  we  summarize  briefly  the  work 
completed  under  the  contract  during  the  year  1  April  1971 
through  31  March  1972.  This  project,  as  the  title  indicates, 
covers  a  wide  range  of  topics.  A  number  of  specific  investi¬ 
gations  were  undertaken  within  the  broad  scope  of  the 
problem  of  discriminating  earthquakes  from  underground 
nuclear  explosions. 

The  topics  studied  can  be  grouped  under  three  general 
headings : 

1.  Wave  propagation  in  laterally  heterogeneous  media, 

2.  Source  mechanisms  of  earthquakes  and  explosions, 

3 .  Earth  structure . 

In  the  following  sections  we  present  abstracts  of  the 
published  papers  falling  in  each  category.  Unpublished 
results  are  discussed  in  more  detail.  Also  included  are 
lists  of  all  publications  and  theses  supported  under  this 
project  during  the  contract  year. 

Much  has  been  learned  about  the  interior  of  the  earth 
by  interpreting  seismic  data  in  terms  of  laterally  homo¬ 
geneous  earth  models.  To  unravel  the  fine  structure  of  the 
crust  and  upper  mantle,  and  to  understand  the  effects  of 
such  fine  structure  on  seismic  waves,  more  realistic,  laterally 
heterogeneous  earth  models  must  be  studied.  Several  such 
investigations  are  summarized  in  this  report. 
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In  his  Ph.D.  thesis  work,  R.  Ward  developed  a  theoretical 
technique  to  compute  the  P  waveforms  from  sources  located 
near  inhomogeneities.  Since  many  inhomogeneous  zones  in  the 
earth  have  characteristic  widths  similar  to  the  wavelength 
of  long-period  body  waves,  this  treatment  allows  wave  ampli¬ 
tudes  to  be  calculated  with  greater  accuracy  than  does  the 
classical  ray-theory  approach.  Specific  problems  to  which 
the  technique  has  been  applied  include  the  synthesis  of  P 
waves  from  sources  near  a  phase-change  boundary  in  the  mantle 
and  near  a  descending  slab  of  lithosphere  as  at  island  arcs. 

Also  as  a  Ph.D.  thesis,  R.  Madariaga  treated  for  the 
first  time  the  theoretical  effect  of  lateral  heterogeneities 
in  the  earth  on  free  oscillations  and  long-period  surface 
waves.  Using  a  perturbation  theory  method,  he  demonstrated 
that  for  realistic  earth  models  with  structurally  distinct 
oceanic  and  continental  regions  there  is  significant 
splitting  of  earth  normal  mode  spectral  peaks  for  angular 
order  numbers  above  a  critical  value,  in  agreement  with 
observations.  Further,  standard  techniques  for  "regionalizing" 
the  earth  from  the  dispersion  of  long-period  surface  waves 
were  brought  into  question. 

Several  approaches  to  the  questions  of  the  forces 
triggering  earthquakes  and  the  physics  of  the  earthquake  source 
are  currently  being  followed.  Recent  data  has  required  some 
revision  of  Aki's  scaling  law  of  far-field  seismic  spectra. 
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The  law  is  still  valid  for  large  magnitudes  (Ms  >  6)  and 
periods  (T  >  10) ,  but  one  must  conclude  that  large  earth¬ 
quakes  and  small  ones  are  fundamentally  dissimilar.  At 
least  two  simple  modifications  of  the  scaling  law  can 
explain  the  data. 

An  alternative  view  of  the  earthquake  source  is  that 
of  Ida,  who  models  a  fault  as  a  rupturing  crack.  Crack 
theory  may  be  used  to  obtain  improved  relationships  for 
rupture  velocity  and  the  stress  field  near  the  crack  for 
later  comparison  of  calculations  with  near-field  seismic  data. 

The  question  of  whether  tides  play  an  important  role  in 
triggering  earthquakes  was  reexamined  by  Shlien.  An 
extensive  comparison  of  earthquake  occurrences  with  tidal 
stresses  was  conducted,  but  no  significant  correlation  could 
be  established. 

For  many  underground  nuclear  explosions,  the  surface 
waves  indicate  a  substantial  double  couple  component  at  the 
source,  most  likely  due  to  release  of  tectonic  strain  in  the 
rock  surrounding  the  explosion.  Toksttz  and  Kehrer  have 
demonstrated  that  such  strain  release  depends  on  the  rock 
type  and  pre-stress.  Tectonic  strain  energy  released  as 
surface  waves  can  occasionally  exceed  the  surface  wave  energy 
from  the  explosion  alone.  For  such  explosions,  however,  the 

M_  -  mu  discriminant  is  still  effective. 

S  b 

Many  aspects  of  earth  structure  need  better  definition 
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if  the  effects  of  propagation  path  and  source-region 
peculiarities  are  to  be  successfully  removed  from  seismic 
signals.  Two  of  these  aspects  have  to  date  received 
insufficient  attention:  seismic-wave  attenuation  and  ambient 
stress.  The  former  quantity  is  important  because  of  its 
effect  on  the  absolute  amplitude  and  spectral  shape  of 
seismic  waves,  the  latter  because  of  its  relation  to 
natural  seismicity  and  earthquake  source  parameters. 

A  detailed  investigation  by  Solomon  of  seismic  attenuation, 
or  Q,  beneath  North  America  has  demonstrated  that  Q  is 
frequency  dependent,  in  agreement  with  predictions  of 
attenuation  in  a  partially  melted  upper  mantle.  Together  with 
the  known  dependence  of  Q  on  depth  and  its  variation  laterally, 
this  fact  complicates  the  correction  of  seismic  wave  spectra 
for  attenuation  and  thus  many  discrimination  criteria  as 
well.  In  particular,  the  Mg  -  it^  discriminant,  the  estimation 
of  source  depth  from  surface  wave  spectra,  and  the  deter¬ 
mination  of  long-period  surface  wave  magnitude  and  other 
source  parameters  can  be  hampered  by  imprecise  or  incorrect 
knowledge  of  Q. 

As  a  step  toward  linking  dynamical  models  of  the 
earth's  crust  and  mantle  with  global  patterns  of  earthquake 
occurrence  and  source  mechanisms ,  Smith  and  Toksfiz  have  calcu¬ 
lated  the  stress  fields  within  geophysically  reasonable 
models  of  descending  lithosphere  slabs.  It  is  within  and 
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adjacent  to  these  slabs  that  a  majority  of  the  world's 
earthquake  occur.  The  models  of  stress  patterns,  including 
the  effects  of  temperature  and  phase  changes,  are  in  good 
agreement  with  the  stress  distributions  inferred  from  the 
focal  mechanisms  of  deep  earthquakes. 
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2.  WAVE  PROPAGATION  IN  LATERALLY  HETEROGENEOUS  MEDIA 

2 . 1  Synthesis  of  Teleseismic  P-Waves  from  Sources  Near 

Transition  Zones  by  Ronald  W.  Ward  (Abstract) 

We  have  developed  a  useful  method  to  compute  both  short- 
period  and  long-period  theoretical  seismograms  from  a  P-wave 
source  in  or  near  an  inhomogeneous  transition  zone.  The 
method  is  applicable  to  the  cases  in  which  the  characteristic 
dimension  of  the  transition  zone  is  comparable  to  the  wave¬ 
length  and  the  source  is  located  near  (within  a  few  wave¬ 
lengths)  or  in  the  transition  zone.  Recent  evidence 
indicates  that  such  transition  zones  may  exist  in  the  upper 
mantle.  The  transition  zone  distorts  the  wave  form,  gives 
rise  to  focusing  and  complex  interference  with  reflected 
waves,  and  also  forms  frequency  dependent  shadow  zones.  The 
method  can  take  a  fuller  account  of  these  effects  than  ray- 
theoretical  methods  do,  in  order  that  the  nature  of  these 
transition  zones  may  be  determined  more  precisely. 

In  our  method  we  avoid  the  complexity  of  scattering  in 
the  nearfield  of  a  point  source  by  considering  the  reciprocal 
problem,  which  is  simply  related  to  the  forward  problem.  We 
solve  the  reciprocal  problem  by  using  a  plane  wave  approximation 
in  the  far  field.  An  exact  wave-solution  for  plane  waves 
incident  on  a  one-dimensional  acoustic  or  elastic  linear 
transition  layer  is  used.  An  algorithm  is  developed  to 
evaluate  the  wave-solutions  for  a  general  class  of  linear 
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transition  zones.  The  frequency  domain  far  field  response 
to  a  source  near  a  transition  zone  is  synthesized  to  obtain 
theoretical  seismograms  for  comparison  with  observations. 

We  apply  this  method  to  the  problem  of  the  20°  dis¬ 
continuity,  by  synthesizing  short-period  and  long-period 
seismograms  for  several  models  of  the  transition  zone.  We 
find  that  the  total  velocity  increase  across  the  trcinsition 
zone  may  be  accurately  determined  from  several  independent 
observable  features  of  the  seismograms.  Even  if  the  focal 
depth  is  not  accurately  known,  one  can  determine  the  thickness 
of  the  transition  zone  from  the  reflections,  provided  there 
is  a  low  noise  (including  corrupting  arrivals)  level  (S/N 
40  db) .  If  the  total  velocity  increase  and  the  focal  depth 
are  precisely  known,  one  can  infer  the  velocity  gradient  in 
the  transition  zone  from  the  observations  of  a  source  within 
the  transition  zone. 

We  use  our  method  to  treat  a  simplified  model  of  the 
sinking  lithospheric  slab  with  a  10%  maximum  velocity  contrast. 
For  this  model,  we  find  that  the  ray-theoretical  solution 
predicts  the  amplitude  of  the  short-period  seismogram  quite 
well.  The  amplitude  of  the  long-period  seismogram,  however, 
does  not  follow  the  ray-theoretical  prediction  well.  For 
example,  the  amplitude  in  the  ray-theoretical  shadow  zone 
may  be  typically  50%  of  the  direct  P-wave  amplitude.  The 
lithospheric  slab  broadens  the  waveform  of  the  long-period 
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P-wave  in  certain  radiation  azimuths.  Source  mechanism 
investigators  working  in  the  time  domain  attribute  the  cause 
of  similar  broadening  of  the  waveform  to  the  finiteness 
of  the  earthquake  source.  We  show  that  for  sources  occurring 
in  and  near  the  lithospheric  slab  the  broadening  of  the 
waveform  by  the  slab  must  be  taken  into  account  in  source 
dynamic  studies. 

We  find  that  the  major  cause  of  the  difference  in  the 
observed  radiation  pattern  between  Longshot  and  Milrow  is 
more  likely  due  to  source  effects  than  due  to  the  effect  of 
a  nearby  lithospheric  slab  on  the  different  dominant  fre¬ 
quencies.  For  example,  the  cause  of  this  amplitude  difference 
may  be  due  to  greater  strain  energy  release  for  Milrow  than 
Longshot,  introducing  a  double  couple  component  into  the 
radiation  pattern. 

2.2  Free  Oscillations  of  the  Laterally  Heterogeneous 

Earth  by  Raul  I.  Madariaga  (Abstract) 

The  splitting  of  the  eigenfrequencies  of  the  earth  due 
to  lateral  heterogeneities  was  studied  by  the  use  of  classical 
theory  of  perturbations  of  the  degenerate  eigenfrequencies 
of  the  spherically  symmetric  average  earth.  The  eigen¬ 
functions  30m(r)  associated  with  a  degenerate  frequency 
n  *< 

lead  assymptotically  (for  large  1)  to  surface  waves  propa¬ 
gating  along  great  circle  paths  determined  by  l  and  m.  The 
lateral  heterogeneities  couple  these  eigenfunctions  and 


produce  a  complicated  interference  pattern  on  the  surface  of 
the  earth.  The  perturbation  expansion  leads,  to  first  order, 
to  a  matrix  eigenvalue  problem  for  a  Hermitian  matrix  which 
selects  the  constructive  interference  patterns.  The  elements 
of  the  matrix  express  the  coupling  between  different  eigen¬ 
vectors  of  the  degenerate  problem. 

We  expand  the  lateral  heterogeneities  in  a  series  of 

spherical  harmonics  <6p  <r)  y  (6  ,ip)  ,  6y  t(r)y  t(d,ip), 

s  s  s  s 

(r)ys  (0,ip)).  Every  term  in  this  expansion  couples  the 
eigenfunctions  selectively.  These  selection  rales  can  all 
be  deduced  from  a  simple  vector  diagram.  If  there  is  no 
strong  coupling  with  overtones  or  spheroidal  -  toroidal 
coupling,  the  odd  terms  ir.  the  spherical  harmonics  expansion 
of  the  heterogeneities  do  not  affect  the  eigenfrequencies. 

For  instance,  there  would  be  no  splitting  if  the  earth  had 
a  continental  hemisphere  and  an  oceanic  hemisphere.  Also, 
terms  such  that  s  >  2 £  do  not  affect  the  eigenfrequencies  of 
modes  of  order  £. 

Observations  indicate  that  splitting  due  to  lateral  hetero¬ 
geneity  is  significant  for  angular  orders  %  >  20  or  T  <  350 
sec  and  is  of  the  order  of  l.%  -  1.5%  for  toroidal  oscillations 
at  H  =  40  and  about  1%  for  spheroidals.  Numerical  computations 
for  proposed  models  of  continental  and  oceanic  structures 
show  that  the  splitting  predicted  theoretically  is  of  the 
order  of  the  observed.  A  model  in  which  lateral  heterogeneity 
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is  caused  solely  by  anomalous  structure  of  the  tectonically 
active  regions  was  also  studied.  Computations  based  on 
realistic  structures  beneath  these  regions  show  that  the 
effect  is  too  small  to  account  for  the  observed  amount  of 
splitting. 

The  great  circle  ray  theory  commonly  used  to  study  the 
dispersion  of  long  period  surface  waves  is  shown  to  be  valid 
only  if  the  heterogeneities  can  be  expanded  with  spherical 
harmonics  of  very  low  s(s  <<  £,)  .  This  puts  a  severe  con¬ 
straint  on  the  applicability  of  this  theory,  in  particular 
it  cannot  be  applied  on  tectonic  regions.  Scattering  and 
interference  of  surface  waves  in  the  laterally  heterogeneous 
earth  produce  time  varying  amplitude  and  phase,  which  cannot 
be  interpreted  as  due  to  the  effect  of  local  structure  along 
the  wave  path. 

2 . 3  Spectral  Splitting  of  Toroidal-Free  Oscillations  Due  to 

Lateral  Heterogeneity  of  the  Earth ' s  structure  by 

Raul  Madariaga  and  Keiiti  Aki  (Abstract) 

Contradicting  results  have  been  obtained  by  H.  Kanamori 
(1970)  and  A.M.  Dziewonski  (1970)  on  the  'pure-path'  dis¬ 
persions  of  mantle  Rayle^g1:  waves.  This  discrepancy  may  be 
due  to  the  inadequacy  of  their  methods  of  ray-theory  inter¬ 
pretation,  in  that  the  total  phase  delay  is  assumed  to  be 
the  algebraic  sum  of  delays  over  the  lengths  of  pure  paths 
measured  along  the  great-circle  path.  In  order  to  take  into 


-11- 


account  the  effect  of  interference  by  the  waves  propagating 
in  other  directions,  we  apply  a  perturbation  method  to 
several  laterally  heterogeneous  earth  models  and  calculate 
the  splitting  of  exgenfrequency  of  toroidal  oscillations. 

For  the  order  number  A  =  40,  observation  shows  the  variation 
of  the  spectral  peak  to  be  as  much  as  1.7%  of  the  degenerate 
frequency.  The  continent-ocean  heterogeneity  with  8099  or 
5.08M  as  the  oceanic  model  and  CANSD  as  the  continent  model 
produces  the  multiplet  width  of  about  1%.  The  continent- 
ocean  heterogeneity  with  01  and  Si  as  the  constituent  models 
gives  less  than  0.4%  multiplet  width.  The  tectonic  hetero¬ 
geneity,  in  which  the  tectonically  active  region  is  assumed 
to  be  different  from  the  rest  of  the  world,  also  produced 
too  small  a  multiplet  width  for  any  choice  of  models  of  these 
regions.  The  ray- theory  method  of  interpreting  the  great- 
circle  phase  velocities  has  been  examined  by  the  criterion 
obtained  by  R.  Madariaga  (1972)  .  With  the  modes  having  order 
numbers  H  <  40,  the  ray  theory  is  not  applicable  to  any  of 
the  laterally  heterogeneous  earth  models  discussed  here. 
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3.  SOURCE  MECHANISMS  OF  EARTHQUAKES  AND  EXPLOSIONS 

3  -1  Earthquake  Mechanism  by  Keiiti  Aki  (Abstract) 

Important  progress  has  been  made  in  the  study  of  the 
earthquake  mechanism  during  the  Upper  Mantle  Project,  through 
the  establishment  of  an  appropriate  mathematical  framework 
which  relates  the  observed  seismogram  with  the  slip  motion 
across  a  fault  plane.  Since  an  arbitrary  fault  slip  is 
described  by  a  function  of  time  and  two  space  coordinates, 
a  complete  inversion  of  the  observed  seismogram  is  practically 
impossible.  The  only  practical  inversion  method  is  to 
describe  the  kin  smatics  of  rupture  growth  along  a  fault 
plane  using  a  small  number  of  source  parameters,  and  then  to 
determine  those  parameters  from  the  seismograms.  A  number 
of  kinematic  models  have  been  proposed  to  cover  a  class  of 
earthquakes.  Preliminary  attempts  have  also  been  made  to 
find  dynamic  models  by  solving  the  problem  of  spontaneous 
rupture  propagation  for  given  initial  conditions. 

There  are  a  few  source  parameters,  which  represent  some 
averaged  quantities  over  the  source  time  and  source  space 
and  are  therefore  indepen^  ent  o?  the  details  of  rupture 
kinematics.  They  are:  (1)  the  seismic  moment,  which  is 
proportional  to  the  total  displacement  averaged  over  the  fault 
surface?  and  (2)  the  apparent  stress,  which  is  proportional 
to  the  ratio  of  total  seismic  energy  to  seismic  moment.  In 
addition  to  these,  the  average  stress  drop  can  be  approximately 
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calculated  from  the  seismic  moment  and  the  fault  area,  without 
detailed  knowledge  of  the  kinematics.  The  accuracies  of 
determining  these  average  quantities  are  discussed  for  the 
Parkfield,  California  earthquake  of  1966,  one  of  the  best 
studied  earthquakes  during  the  U.M.P.  period.  Recent  results 
from  several  well-studied  earthquakes  are  then  summarized  as: 

(1)  the  relation  between  magnitude  and  seismic  moment;  and 

(2)  the  relation  between  stress-drop  and  magnitude.  It  was 
found  that  the  stress  drop  in  shallow  earthquakes  is  10-100 
bar,  apparently  independent  of  magnitude  for  M  >  6.  Finally, 
recent  evidences  for  the  fault  origin  of  deep-focus  earth¬ 
quakes  are  presented. 

3 . 2  Recent  Results  on  the  Mechanism  of  Earthquakes  with 

Implications  on  the  Prediction  and  Control  Program 

by  Keiiti  Aki  (Abstract) 

The  current  earthquake  prediction  and  control  research 
programs  are  based  on  the  working  hypothesis  that  an  earth¬ 
quake  is  a  release  of  tectonic  strain  energy  stored  up  in  a 
volume  of  the  earth's  crust.  Accumulating  seismological 
evidence  indicates  that  this  energy  release  is  done  by  a 
propagating  rupture  which  creates  a  shear  dislocation  over 
a  finite  fault  plane. 

The  velocity  of  the  rupture  propagation  is  the  single, 
most  important  factor  in  the  dynamics  of  energy  release  in 
this  process.  If  the  velocity  is  very  slow  and  the  dislocation 
takes  place  quasi-statically ,  all  the  strain  energy  will  be 


spent  against  friction  on  the  fault  and  there  will  be  no 
seismic  energy  radiation  and  no  earthquake  hazard. 

If  the  velocity  of  rupture  propagation  is  infinite  and 
the  stress  release  takes  place  instantaneously  over  the 
entire  fault  plane,  then  the  particle  velocity  on  the  fault 
is  simply  the  stress  divided  by  the  impedance,  that  is,  the 
product  of  density  and  shear  velocity.  As  discussed  by 
Brune,  the  particle  velocity  of  about  1  meter/sec  is  pre¬ 
dicted  for  the  stress  release  of  about  100  bars. 

For  a  finite  rupture  velocity,  however,  the  impedance 
depends  on  rupture  velocity.  For  example,  for  longitudinal- 
shear  cracks,  the  impedance  becomes  zero  when  the  rupture 
propagates  with  the  shear  velocity.  In  this  case,  the 
particle  motion  becomes  infinite  when  there  is  a  finite 
stress  release  along  the  fault.  For  a  given  tectonic  stress, 
therefore,  one  should  expect  seismic  motion  with  an  amplitude 
from  zero  to  infinity  critically  dependent  on  the  rupture 
velocity.  If  we  find  the  way  to  control  the  rupture  velocity, 
we  can  control  an  earthquake. 

Early  seismological  determination  of  rupture  velocities 
indicated  the  velocity  at  3  to  4  km/sec,  a  little  below  the 
medium  shear  velocity.  Later  studies  extended  the  range  from 
about  1.5  to  about  4 . 5  km/sec . 

We  have  not  yet  established  any  systematic  dependence 
of  rupture  velocity  on  the  earthquake  parameters  such  as 


focal  depth,  fault-plane  geometry,  and  other  tectonic 
features  of  an  epicenter.  For  example,  a  typical  strike-slip 
earthquake  on  the  San  Andreas  fault  and  a  deep  South  American 
earthquake  share  the  same  rupture  velocity  of  2.2  km/sec. 

There  are  some  indications  that  the  rupture  propagation  is 
faster  for  large  earthquakes  than  for  small  ones,  but  we 
do  not  have  enough  data  to  draw  definite  conclusions. 

3.3  Cohesive  Force  across  the  Tip  of  a  Longitudinal-Shear 

Crack  and  Griffith's  Specific  Surface  Energy  by 

Yoshiaki  Ida  (Abstract) 

The  cohesive  force  across  the  fault  plane  is  considered 
in  order  to  understand  the  physical  mechanism  of  rupture  at 
the  tip  of  a  longitudinal-shear  crack.  The  elastic  field 
around  the  tip  of  a  crack  and  the  condition  of  rupture  growth 
are  systematically  derived  from  the  assumption  that  the 
cohesive  force  is  given  as  a  function  of  the  displacement 
discontinuity.  This  assumption  is  more  physically  meaningful 
than  those  originally  used  by  6.1.  Barenblatt  in  1959  and  1962. 
The  stress  field  around  the  tip  is  calculated  for  several 
models  of  cohesive  force,  and  is  shown  to  be  nonsingular  even 
at  the  tip.  The  condition  of  rupture  growth  that  is  used 
to  determine  the  rupture  velocity  turns  out  to  be  equivalent 
to  the  Griffith  criterion  and  the  relation  employed  by  B.V. 
Kostrov  in  1966,  but  the  specific  surface  energy  is  defined 
more  clearly  in  this  paper. 


3.4  Earthquake-Tide  Correlation  by  Seymour  Shlien  (Abstract) 

The  detection  of  moonquakes  that  occur  when  the  moon 
is  at  perigee  has  prompted  a  search  for  tidal  effects  on 
earthquake  occurrences .  An  attempt  was  made  to  correlate 
earthquakes  listed  in  the  CGS-NOAA  epicenter  determinations 
with  the  tidal  phase  of  semidiurnal  tides.  This  study  was 
confined  to  several  seismic  regions  representative  of 
tectonic  and  non-tectonic  regions.  An  extended  form  of 
Schuster's  test  was  used  to  decide  whether  significant 
correlations  existed.  Though  some  tidal  influences  could 
be  accepted  at  a  5%  significance  level,  the  effect  was  not 
consistent  or  stable  with  time.  Earthquakes,  if  they  are 
affected  by  tides,  show  a  slight  tendency  to  occur  at  times 
when  the  tidal  stress  is  changing  most  rapidly.  Insufficient 
data  was  available  to  compare  tectonic  to  non-tectonic  areas. 

An  analysis  of  the  Japanese  aftershock  sequence  which 
began  11  August  1959  was  found  to  have  no  significant  tidal 


correlation. 


3.5  Scaling  Law  of  Earthquake  Source  Time-Function  by 

Keiiti  Aki 

Summary 

Further  evidences  support  the  scaling  law  of  far-field 
seismic  spectrum  based  upon  the  w-square  model  (Aki,  1967) 
for  earthquakes  with  Mg  >  6  and  for  periods  T  >  10  seconds. 
Recent  observations,  however,  unequivocally  require  the 
modification  of  the  above  law  for  periods  T  <  10  seconds. 
Unfortunately,  the  presently  available  data  are  not  sufficient 
for  a  unique  revision  of  the  scaling  law.  We  propose  two 
alternatives  and  discuss  their  implications  and  consequences. 
In  either  case,  we  have  to  conclude  that  a  large  earthquake 
and  a.  small  one  are  substantially  different.  One  interesting 
feature  of  the  w-square  model  appears  to  be  unaffected  by 
the  required  revision;  that  is,  the  spectral  density  of  the 
dislocation  time-function  for  periods  T  <  5  seconds  takes  the 
same  absolute  value,  independent  of  magnitude,  for  earthquakes 
greater  than  M  =6.5.  This  result  has  important  consequences 
in  earthquake  engineering  because  the  seismic  motion  in  the 
vicinity  of  an  earthquake  fault  will  scale  as  the  dislocation 


function. 


1 .  INTRODUCTION 


As  a  first  step  to  relate  the  seismic  spectrum  with 
earthquake  magnitude,  a  model  of  earthquake  ensemble  was 
proposed  by  Aki  (1967)  on  the  basis  of  a  dislocation  theory 
of  earthquake  faulting.  In  this  model,  the  source  factor  of 
far-field  spectrum  diminishes  inversely  proportional  to  the 
square  of  frequency  co  beyond  a  corner  frequency.  For  this 
frequency  dependence,  it  was  called  the  "w-square  model."  Be¬ 
low  the  corner  frequency,  the  spectrum  is  flat  with  the  height 
proportional  to  the  seismic  moment  (Aki  1966) .  A  family  of 
such  spectral  curves  was  constructed  on  the  assumption  that 
large  and  small  earthquakes  are  similar  phenomena  in  a  medium 
with  given  elastic  constants  and  density.  The  assumption  implies 
the  same  geometry,  a  constant  stress-drop,  constant  rupture 
velocity  and  slip  velocity,  independent  of  magnitude,  and  it 
follows  that  the  corner  frequency  is  inversely  proportional 
to  the  fault  length,  and  the  seismic  moment  to  the  cube  of 
fault  length.  Thus,  the  corner  frequency  lies  on  a  straight 
line  with  slope  3  in  Fig.  1  which  shows  the  logarithm  of 
spectrum  against  the  logarithm  of  period.  The  spacings 
between  the  curves  are  made  equal  at  the  period  of  20  seconds 
to  be  compatible  with  the  definition  of  M  by  Gutenberg  and 


Richter.  The  scaling  law  shown  in  Fig.  1  explained  very 
well  Berckhemer's  (1962)  observations  on  spectral  ratios, 
Gutenberg-Richter ' s  (1956)  Ms-Mb  relation  for  Ms  >  6, 
and  Tocher's  (1960)  data  on  the  earthquake  fault  length  and 
magnitude . 

The  assumption  of  similarity  comes  from  the  idea  advo¬ 
cated  by  Tsuboi  (1940,  1956,  1958,  1965)  who  held  the  view 
that  the  strain  energy  density  prior  to  the  earthquake  occur¬ 
rence  is  the  property  of  rock  material  independent  of  earth¬ 
quake  magnitude,  and  that  the  energy  of  an  earthquake  is  de¬ 
termined  by  the  volume  within  which  it  has  been  stored.  This 
idea  radically  contradicts  the  assumption  underlying  Benioff's 
(1951)  strain  release  curve,  which  was  calculated  as  the  square 
root  of  energy  release  implying  a  constant  earthquake  volume 
independent  of  magnitude.  The  latter  view  may  be  natural 
to  the  California  seismologists,  becuase  the  majority  of  Cali¬ 
fornia  earthquakes,  large  or  small,  appear  to  be  associated 
with  the  same  fault  plane;  the  San  Andreas.  It  was,  however, 
unacceptable  to  the  Japanese  seismologists  who  were  familiar 
with  the  distinct  difference  in  spectral  structure  between 
large  and  small,  especially  microearthquakes  (Asada  1957) . 

The  controversy  appeared  to  have  been  settled  when  B&th  and 
Ouda  (1964)  summarized  observational  data  in  favor  of  Tsuboi' s 
idea  and  proposed  an  improvement  of  Benioff's  method  for 


calculating  the  strain  release  curve.  ri?hus,  the  assumption 
of  similarity  was  the  most  reasonable  one  at  the  time  when 
the  w-square  model  was  proposed. 

Since  then,  a  large  amount  of  observational  data  have 
become  available  for  more  critical  testing  of  the  proposed 
model.  In  general,  the  new  evidences  support  the  w-square 
model  for  earthquakes  larger  than  Ms  =  6,  and  for  periods 
longer  than  10  seconds.  For  smaller  earthquakes  or  for 
shorter  oeriods,  new  evidences  require  revision  of  the  w-square 
model.  The  purpose  of  the  present  paper  is  to  propose  such  a 
revision,  and  to  discuss  its  implications  and  consequences. 
Since,  unfortunately,  the  presently  available  data  are  not 
sufficient  for  a  unique  revision,  discussion  will  be  made  for 
two  alternative  revision. 

An  interesting  feature  of  the  u)-square  model  is  that  the 
spectral  density  of  dislocation  time-function  for  periods 
shorter  than  5  seconds  is  identical,  in  the  absolute  scale, 
for  all  the  earthquakes  greater  than  Mg  =  6.5.  This  feature 
of  the  w-square  model,  which  is  particularly  important  for 
earthquake  engineering,  is  unaffected  by  the  revision  required 
in  the  present  paper. 
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2.  COMPARISON  WITH  OBSERVATIONS 

The  validity  of  the  scaling  law  shown  in  Fig.  1  can  be 
tested  against  various  kinds  of  observations.  In  general, 
the  right-half  (T  >  10  sec)  of  Fig.  1  shows  excellent  agree¬ 
ments  with  the  observed  data,  but  the  left-half  (T  <  10  sec) 
does  not.  Let  us  start  vrith  observations  at  the  infinite 
period,  that  is,  static  fault  parameters. 

(1)  Fault  length .  Fig.  2  shows  the  data  on  magnitude 
and  fault  length  reproduced  from  Chinnery  (1969) .  In  view 
of  the  difficulties  in  extracting  reliable  observations  of 
fault  parameters  from  field  data  the  observed  points  in  Fig. 

2  were  thoughtfully  limited  to  those  of  near-vertical  faults 
on  which  the  movement  was  predominantly  strike-slip.  The 
surface  wave  magnitude  Mg  is  used  for  events  with  magnitude 
over  6,  and  local  Richter  magnitude  for  events  less  than 

6.  The  empirical  formulas  by  Iida  (1959,  1965)  and  Tocher 
(1958)  show  a  good  fit  to  the  data  for  large  earthquakes. 
Otsuka's  (1965)  formula  is  proposed  to  take  into  account  the 
hidden  part  of  a  fault  unnoticed  by  field  observation.  Both 
Otsuka's  and  Press'  (1967)  formula  predicts  much  steeper  slope 
than  observed  for  V  ^  7.  This  discrepancy  was  attributed  to 
the  non-similar  fault  shape  between  large  and  small  earthquakes 
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due  to  the  effect  of  crustal  structure. 

The  a)- square  model ,  on  the  other  hand,  gives  exactly 
the  same  slopes  as  those  of  the  Iida  and  Tocher  formulas,  and 
explains  the  observed  slope  for  large  events  excellently  with¬ 
out  invoking  the  non-similarity.  The  bending  of  the  curve 
at  about  M  =  7  is  due  to  the  inefficiency  of  Mg  as  a 
measure  of  the  size  of  large  earthquakes. 

The  curve  for  the  w-square  model  in  Fig.  2  was  drawn, 
first,  by  finding  the  relation  between  Ms  and  the  corner 
period  from  Fig.  1,  and  then  multiplying  the  corner  period 
by  a  constant  to  obtain  the  fault  length.  Under  the  assump¬ 
tion  of  similarity,  the  corner  period  and  fault  length  should 
be  proportional  to  each  other.  The  proportionality  constant 
giving  the  best  fit  to  observation  is  0.65  km/ sec  (see  Fig. 
11  of  Aki  1967)  .  It  is  interesting  to  compare  this  value  with 
the  theoretical  coefficient  used  by  Brune  (1970,  1971)  ,  in 

which  the  radius  r  of  earthquake  source  is  related  to  the 

2.343 

corner  frequency  fo  of  shear  wave  spectra  by  r  =  2-nfQ  ' 
where  3  is  the  shear  velocity.  If  we  take  2r  as  the  fault 
length,  the  coefficient  is  2.6  km/sec  for  3  =  3.5  km/sec. 
This  is  about  4  times  larger  than  the  value  obtained  from  the 
observational  data  for  large  earthquakes.  It  appears  that 
Brune' s  theory  does  not  apply  to  large  earthquakes.  The  dis- 
be  attributed  to  his  assumption  of  infinite 


crepancy  may 
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rupture  velocity,  which  is  not  realistic  for  a  spontaneous 
rupture  (Ida  &  Aki  1972) . 

Within  the  scheme  of  Aki's  (1967)  statistical  fault 
model,  the  corner  period  T  is  related  to  the  mean  free  path 
kL'^  of  fault  propagation  by  the  relation 


where  v  is  the  velocity  of  rupture  propagation.  For  v  =  3  km, 
we  find  that  the  mean  free  path  is  proportional  to  the  corner 
period  with  coefficient  about  0.5,  which  is  close  to  the  observed 
coefficient  0.65  for  the  fault  length.  For  large  earthquakes 
and  for  long  periods,  therefore,  the  rupture  propagation  appears 
to  be  smooth  and  encounters  no  obstacle  during  the  growth  to 
its  final  length. 

The  curve  based  on  the  w-square  model  does  not  explain  the 
data  for  small  earthquakes  shown  in  Fig.  2,  which  includes  the 
Imperial  earthquake  of  March  4,  1966  described  by  Brune  and 
Allen  (1967)  who  demonstrated  beyond  doubt  that  an  earthquake 
with  magnitude  3.6  could  have  a  10  km  fault  length.  However, 
some  of  the  data  are  questionable.  For  example,  one  after¬ 
shock  (M  -  4.9)  of  the  Parkfield  earthquake  of  1966  was  given 
a  fault  length  of  33  km  from  creep  observations  (Wyss  &  Brune 
1968) .  The  wide-band  spectra  of  Love  waves  (Filson  &  McEvilly 
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1967)  from  this  earthquake,  however,  did  not  show  the  evi¬ 
dence  for  such  a  long  fault.  Since  the  field  measurement  of 
fault  length  becomes  increasingly  difficult  with  decreasing 
magnitude,  we  must  resort  to  indirect  methods.  For  example, 
Lieberman  and  Pomeroy (1970)  used  the  data  on  an  aftershock 
area  in  a  general  support  of  the  Wyss-Brune  curve  of  Fig.  2. 
However,  an  accurate  determination  of  aftershock  area  for  a 
small  main-shock  is  a  difficult  problem  and  also  the  after¬ 
shock  area  may  give  an  overestimate  of  the  main-shock  fault 
area  (Aki  1968) . 

(2)  Stress  drop.  Further  evidences  from  field  data  support 
the  validity  of  the  w-square  model  for  large  earthquakes. 

Fig.  3  shows  the  relation  between  Ms  and  log  LD2  reproduced 
from  King  and  Knopoff  (1968) ,  where  L  is  the  fault  length 
and  D  is  the  fault  offset.  King  and  Knopoff  found  that 
the  slope  of  best-fitting  lines  is  significantly  different 
from  the  slope  of  any  published  magnitude-seismic  energy  re¬ 
lations.  From  this  result,  assuming  a  dislocation  model  with 
a  conrtant  efficiency  independent  of  magnitude,  they  concluded 
low  fractional  stress-drop  for  small  earthquakes.  Their  con¬ 
clusion  has  been  extended  to  small  earthquakes  by  Wyss  (1970) , 
and  was  considered  as  an  evidence  against  the  similarity  assump¬ 
tion  (constant  stress-drop)  underlying  the  w-square  model. 
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2  3 

The  parameter  LD  is  proportional  to  L  ,  and  there¬ 
fore  to  the  seismic  moment  in  the  u-square  model.  One  can 

2 

draw  a  theoretical  relation  between  Mg  and  LD  by  finding 

the  seismic  moment  for  a  given  M  in  Fig.  1  and  multiplying 

s 

a  constant.  For  a  rectangular  fault  with  width  W,  this 

constant  is  equal  to  which  is  a  measure  of  strai.i- 

^  -17 

release  or  stress  drop.  The  constant  was  chosen  as  3  x  10 

11  —2 

dyne-1  cm2  in  Fig.  3.  Assuming  that  p  =  3  x  10  dyne  cm 
in  the  earth's  crust,  we  find  that  ^  =  10  5  and  =  30  bars. 
This  value  is  a  reasonable  rough  estimate  of  stress  drop  in 
large  shallow  earthquakes  (Chinnery  1964;  Brune  &  Allen  1967). 

A  recent  summary  of  earthquake  mechanism  studies  by  Aki  (1972) 
demonstrates  that  the  stress  drop  in  shallow  earthquakes  with 
M  >  6  is  10  to  100  bars  independent  of  magnitude.  Here,  again, 
we  see  that  the  oi-square  model  can  explain  observed  field  data 
without  invoking  non-similarity . 

(3)  Seismic  moment.  Fig.  4  is  reproduced  from  Aki  (1972) , 

who  summarized  the  data  on  M  and  seismic  moment  obtained 

s 

from  long-period  surface  waves  and  free  oscillations.  Body 
wave  results  were  not  included  because  of  the  controversial 
window  effect  discussed  by  Linde  &  Sacks  (1971) ,  who  con¬ 
cluded  that  theories  which  predict  constant  displacement 
spectrum  for  body  waves  at  long  periods  (the  w- square  model 
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is  one  of  them)  must  be  in  error.  Their  results  are  in  favor 
of  the  Archambeau  (1968)  theory  which  predicts  a  sharp  drop  of 
spectrum  toward  zero  frequency.  We  disagree  with  their  con¬ 
clusion  because  the  surface  wave  and  free  oscillation  spectra 
which  do  not  suffer  from  the  window  effect  invariably  show 
consistency  with  the  assumption  of  step  function  dislocation 
at  the  low-frequency  end  of  observable  spectrum.  The  best 
example  is  the  free  oscillation  amplitude  excited  by  the  Kurile 
earthquake  of  October  13,  1963  (M  =  81/4).  Abe  (1970)  demon¬ 
strated  that  the  assumption  of  step  function  was  valid  for  the 
order  numbers  10  to  30  for  both  spheroidal  and  torsional  oscil¬ 
lations.  Other  convincing  cases  may  be  found  in  Ben-Menahem 
and  ToksOz  (1963),  Aki  (1966),  Kanamori  (1970a, b) ,  Tsai  and  Aki 
(1969,  1970a, b,  1971)  . 

2 

The  smooth  curve  designated  as  "ui  -model"  in  Fig.  4  was  drawn  by 
plotting  the  height  of  the  flat  portion  of  spectrum  for  a 
given  Mg  shown  in  Fig.  1.  The  absolute  value  was  fixed  in 
such  a  way  that  the  curve  passes  through  the  observed  point  for 
the  Niigata  earthquake,  for  which  the  first  accurate  determin¬ 
ation  of  seismic  moment  was  made  by  Aki  (1967)  . 

The  dashed  lines  in  Fig.  4  represent  the  calibration  curve 
used  by  Brune  (1968)  and  Davies  and  Brune  (1971)  to  find  the 
rate  of  slip  along  major  fault  zones  of  the  earth  from  the 
earthquake  magnitude  data.  The  curve  was  developed  through  a 


-27- 


series  of  papers  by  Brune  and  his  colleagues  (Brune  &  King 
1967;  Brune  &  Engen  1969;  Wyss  &  Brune  1968),  and  consisting 
of  the  following  line  segments : 

(A)  log  Mq  =  Ms  +  19.9  for  Ms  >  7.5  and  strike  slip. 

(0.3  is  added  for  dip  slip) 

(B)  log  Mq  =  Ms  +  19.2  for  7  >  Ms  >  6 

(C)  log  Mq  =  1.4Ml  +  17.0  for  6  >  ML  >  3 

The  segments (A) and (B) were  originally  proposed  by  Brune  and 

King  as  a  relation  between  Ms  and  the  Rayleigh  wave  amplitude 

at  about  100  second  periods.  They  assumed  aT1  dependence 
-2 

instead  of  <u  ,  and  using  Tocher  and  Xida's  data  on  fault 
length,  assumed  that  the  corner  period  is  longer  than  100 
seconds  for  Ms  >  7.5,  and  is  shorter  than  20  seconds  for 
Ms  <  7.  Thus,  the  constants  in  (A)  and  (B)  differ  by  the 
logarithm  of  the  ratio  of  two  periods,  as  can  be  figured  out 
from  a  schematic  illustration  of  scaling  law  for  the  "a)-model" 
shown  in  Fig.  5.  It  is  obvious  from  Fig.  5  that  the  coefficient 
in  the  moment-magnitude  relation  for  M  >  7  should  be  1.5 
and  cannot  be  unity  as  given  in  (A) .  In  order  to  be  consistent 
with  the  model  described  by  Brune  and  King  (1967) ,  we  must 
replace  the  formula  (A)  by  the  following : 


tf' 


-28- 


(A1)  log  MQ  =  1.5  Ms  +  16.0,  for  M  >  7 

This  will  connect  to  the  line  (B)  properly  as  expected  for 
the  w-model  sketched  in  Fig.  5. 

Since  the  formula  (C)  was  obtained  empirically  as  the 

relation  for  local  magnitude  Ml,  it  is  not  applicable  to  the 

data  shown  in  Fig.  4.  The  formula  (B)  should  apply  to  all 

Ms  less  than  6  as  can  be  seen  from  Fig.  5.  If  we  compare 

the  data  in  Fig.  4  with  the  formulas  (A*)  and  (B)  shown  by 

chained  lines,  we  find  a  systematic  discrepancy;  all  the  data 

for  M  <6.5  fall  below  the  line  (B) ,  and  all  the  data  for 
s 

Ms  >  7.5  lie  above  the  line  (A*),  except  for  the  Sanriku 
earthquake,  which  is  an  unusual  "lithospheric  normal  faulting 
according  to  Kanamori  (1971)  .  On  the  other  h?.nd,  the  w2-model 
explains  the  observation  except  for  the  Sanriku  earthquake 
without  significant  systematic  error. 

It  is  true,  however,  that  Brune-King's  w-model  explains 
the  mantle  wave  data  of  Brune  and  King  (1967)  and  Brune  and 
Engen  (1969)  somewhat  better  than  the  w-square  model.  Most 
of  their  data  on  large  earthquakes  are,  however,  from  old 
instruments,  and  Brune  and  Engen  express  some  concern  about 
uncertain  instrument  calibration.  Furthermore,  the  correction 
for  attenuation  and  geometrical  spreading  may  be  biased  because 
most  of  the  data  for  small  earthquakes  are  from  the  Pasadena 
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station  and  larger  ones  are  from  other  stations.  Such 
problems  may  be  avoided  by  the  use  of  spectral  ratio  between 
earthquakes  with  different  size  but  with  common  path  and 
recording  station. 

(4)  Spectral  ratio.  The  observed  seismic  spectrum  is  a 
function  of  source,  path  and  receiver.  The  simplest  way  of 
isolating  the  source  spectrum  is  to  compare  seismograms  obtained 
by  the  same  seismograph  at  the  same  station  from  two  earthquakes 
of  the  same  epicenter  but  of  different  size.  Berckhemer  (1962) 
was  able  to  collect  6  such  earthquake  pairs  from  the  Stuttgart 
records  for  the  period  1931  to  1951.  Their  magnitudes  range 
from  4.5  to  8.0,  and  the  epicentral  distances  400  to  9000  km. 

The  separation  between  epicenters  of  each  pair  was  less  than 
1°.  Their  locations  are  shown  in  a  world  map  in  Pig.  6.  The 
w-square  model  explains  very  well  the  observed  spectral  ratio 
for  all  the  pairs,  except  the  one  in  the  Alps  for  which  the 
comparison  was  fair  (Aki  1967) .  A  good  agreement  was  obtained 
also  for  two  aftershocks  of  the  Kern  County,  California  earth¬ 
quake  of  1952. 

A  similar  collection  of  spectra  for  earthquake  pairs  has 
been  made  recently  by  Tsujiura  (1972) ,  using  his  multi¬ 
channel  band-pass  seismographs  (Tsujiura  1966,  1967,  1969). 

The  location  of  earthquake  pairs  is  shown  in  Fig.  6.  The 
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w-square  model  explained  the  spectral  ratio  for  the  pairs  in 
China  (M7.5/M6.1)  and  Halmahera  (M7.2/M6.3)  excellently. 

However,  a  striking  discrepancy  was  found  for  two  pairs  of 
earthquakes  in  the  Aleutians  (M7.1/M6.3,  M7.0/M5.7).  Judging 
from  their  epicenters,  they  all  belong  to  the  underthrusting 
rather  than  the  extensional  group  (Stauder  1968).  Tsujiura's 
observation  indisputably  shows  that  the  M7  and  M6  Aleutian 
earthquake  share  nearly  identical  spectral  shape  in  the 
period  range  10  to  100  seconds.  This  result  definitely  con¬ 
tradicts  the  prediction  of  the  w-square  model,  because,  as 
shown  in  Fig.  1,  the  corner  period  is  about  10  seconds  for 
M=6,  and  about  40  seconds  for  M=7 . 

So  far,  the  failure  of  the  w-square  model  in  explaining 
the  scale  effect  on  seismic  spectrum  may  be  considered  as  an 
exception.  Such  exceptional  cases  are  the  seismic  moment  of 
the  Sanriku  earthquake  and  the  spectral  ratios  of  the  two 
earthquake  pairs  in  the  Aleutians.  So  far,  we  have  considered 
the  period  range  longer  than  10  seconds.  Once  we  enter  the 
period  range  shorter  than  10  seconds,  however,  the  failure 
of  the  w-square  model  becomes  a  rule  rather  than  an  exception. 

(5)  Ms  -  Mb  relation.  Gutenberg  and  Richter  (1956)  discovered 
a  discrepancy  between  the  magnitude  scale  based  upon  short-period 
body  waves  (M,  )  and  that  based  upon  long -period  surface  waves  (Mg) . 
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This  discrepancy  was  attributed  to  the  scale  effect  by  Aki 
(1967)  ,  and  it  was  shown  that  the  w-square  model  explains 
the  Ms  -  Mb  relation  observed  by  them. 

Recently ,  the  Ms  -  Mb  relation  attracted  the  attention 
of  several  seismologists  because  of  its  power  as  a  discrimi¬ 
nant  between  earthquakes  and  underground  explosions  (Lieberman 
et  al.  1966;  Lieberman  &  Pomeroy  1967,  1969;  Capon  et  al.  1967; 
Basham  1969;  Molnar  et  al.  1969;  Evernden  et  al.  1971;  Lieberman 
&  Basham  1971) . 

Fig.  7  summarizes  the  range  of  observed  data  on  Ms  versus 
for  various  regions.  The  curve  designated  as  "w-square 
model"  is  drawn  assuming  that  Ms  is  proportional  to  the  spectral 
density  at  the  period  20  seconds,  and  M5  is  proportional  to 
that  at  1  second.  The  former  assumption  is  valid  because  of 
the  definition  of  Ms,  and  the  latter  is  valid  because  the 
repsonse  of  seismographs  used  for  teleseismic  P  waves  from 
small  events  are  usually  sharply  peaked  at  about  1  second.  This 
assumption  was  used  by  Lieberman  &  Pomeroy  (1969)  in  their 
discussion  of  the  Mg^  relation.  Aki  (1967)  took  into  account 
the  small  effect  of  signal  duration  as  shown  in  Fig.  8  for 
larger  events,  but  such  correction  may  be  unnecessary  for 
smaller  events ,  for  which  the  duration  is  probably  determined 
by  the  instrument  response  and  path  effect.  The  curve  for  the 
w-square  model  is  further  restricted  to  pass  through  the  point 


X 
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(MS  =  =  6  3/4)  according  to  the  original  definition  by 

Gutenberg  &  Richter  (1956) . 

Fig.  7  clearly  demonstrates  that  the  observed  Mg-M^ 
relation  deviates  systematically  from  the  predicted  for  the 
w-square  model.  The  data  follows  the  straight  linear  extra¬ 
polation  of  the  Gutenberg-Richter  formula  (Mg  =  1.59  _  3.97) 

2 

rather  than  following  the  bended  curve  of  the  w  -model  as 
shown  in  Fig.  8.  It  is  extremely  interesting  to  note  that 
this  departure  from  the  w-square  model  makes  the  discrimination 
between  explosion  and  earthquake  possible  for  small  events. 

The  data  for  explosions  are  represented  by  a  line  Mg  =  -  1.8, 

which  was  established  by  Thirlaway  &  Carpenter  (1966) . 

The  fact  that  the  Mg-M^  relation  for  earthquakes  does  not 
bend  according  to  the  w-square  model  but  follows  the  straight 
line  along  the  extrapolated  Gutenberg-Richter  formula  has  an 
additional  support  from  a  work  on  spectral  densities  done  in 
the  U.S.S.R.  Chalturin  (1970)  made  observations  at  Garm, 
using  a  multi-channel  band-pass  seismograph  similar  to 
Tsujiura's  (1966,  1967,  1969).,  and  found  that  the  observed  re¬ 
lation  between  the  amplitudes  at  1  second  and  20  seconds  follows 
a  straight  line  with  the  coeeficient  the  same  as  in  the  Gutenberg 
Richter  formula,  and  does  not  agree  with  the  prediction  of  the 
w-square  model.  This  evidence  is  particularly  strong  because 
no  assumptions  are  made  on  the  relation  between  magnitudes  and 
spectral  densities. 

Evernden  etal.  (1971)  emphasize  the  parallelism  of  the 
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Mg  versus  M^  curves  for  explosions  and  earthquakes,  and 
suggest  that  the  Ms  “  ^  relation  for  earthquakes  with 

<  5  1/4  has  a  slope  of  approximately  1.  Such  a  suggestion 
is  not  inconsistent  with  the  observation,  but  the  data  show  too 
much  scatter  to  allow  a  firm  conclusion  (see  Figs.  7  and  8) . 


3 .  REVISED  MODELS 


Since  the  co-square  model  explains,  in  general,  the  observa¬ 
tions  for  the  period  range  longer  than  10  seconds,  the  right 
half  of  Fig.  1  should  be  left  unchanged.  In  order  to  revise 

the  left  half,  we  shall  consider  the  following  two  extreme 

_2 

cases.  In  one  revision,  we  shall  keep  the  w  -dependence,  but 

discard  the  similarity  assumption  and  change  the  relation  between 

M  and  the  corner  frequency.  This  shall  be  called  the  revised 

-2 

model  A.  In  the  other,  we  shall  give  up  the  w  -dependence  and 
adopt  the  oj-1 -dependence  in  the  period  range  between  10  to  .01 
seconds.  In  this  case,  the  relation  between  M  and  the  corner 
period  is  unchanged  from  the  w-square  model.  We  shall  refer  to 
this  as  the  revised  model  B. 

(1)  Model  A 

Fig.  9  shows  the  family  of  spectral  curves  for  the  revised 

model  A,  in  which  the  spacing  between  the  curves  of  Fig.  1 

at  1  second  was  widened  to  satisfy  observed  Mg  -  M^  relation 

-2 

without  changing  the  w  -dependence.  This  resulted  in 


a  shift  of  corner  period  to  the  right  for  Mg  <  6 . 

For  M  <  4  1/2,  the  corner  period  stays  constant,  because 
s 

we  follow  the  suggestion  of  Evernden  et  al.  (1971)  that  the 
slope  of  the  Mg  -  ^  relation  is  unity  for  ^  <  5  1/4, 

The  relation  between  M  and  the  corner  frequency  is 

5 

shown  in  Fig.  10  for  various  cases.  The  one  corresponding 
to  the  suggestion  of  Evernden  et  al.  indicates  the  corner 
period  of  about  6  seconds  for  Mg  <  4  1/2.  The  physical 
picture  behind  this  is  very  simple.  The  fault  length, 
rupture  propagation  time,  and  rise  time  are  common  to  all 
the  earthquakes  smaller  than  Mg  <  4  1/2.  The  only  difference 
between  them  is  the  offset,  or  strain  release  or  stress  drop. 
This  is  a  revival  of  Benioff's  (1951)  idea  mentioned  in  the 
introduction,  and  was  suggested  by  Brune  &  Wyss  at  the  Woods 
Hole  Conference  on  Seismic  Discrimination,  July,  1970, 
sponsored  by  the  Advanced  Research  Projects  Agency.  This 
idea  is  quite  compatible  with  the  concept  of  "pre-existing 
fault”  which  naturally  violates  the  assumption  of  similarity. 

One  cannot,  however,  extend  the  6  second  corner  period 
to  microearthquakes  (M  =  0) .  The  frequency  ranges  of  usual 
microearthquakes  are  between  10  to  100  cps  [Nevada  may  be  an 
exceptional  area  according  to  Douglas  &  Ryall  (  1972) ,  and 
Takano  (personal  communication)].  For  example,  typical 
records  of  micrearthquakes  may  be  found  in  Furuya  (1969) , 
and  an  extensive  collection  of  data  in  Terashima  (1968)  . 


i 
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A1 though  Mg  cannot  be  used  for  microearthquakes,  the 

seismic  moment  has  been  estimated  for  the  magnitude  zero 

earthquake  by  Wyss  &  Brune  (1968)  ,  Aki  (1969) ,  Takano  (1971) 

and  Douglas  &  Ryall  (1972)  as  1016  ~  1017  dyne  cm.  Therefore, 

the  actual  corner  period  should  reach  the  lower-left  corner 

in  Fig.  10.  The  curve  for  the  w-square  model  does,  but  others 

have  to  be  sharply  bent  to  reach  that  corner.  Such  a  sharp 

bend  in  the  corner  period  curve  produces  a  strange  scaling  law 

of  spectrum  for  Mg  <  3  as  shown  in  Fig.  9.  All  the  earthquakes 

with  -1  <  M  <3  share  the  same  spectral  density  for  frequencies 
s 

higher  than  20  cps. 

Both  the  bend  in  the  magnitude-corner  period  relation  and 
the  above  mentioned  peculiar  scaling  law  of  microearthquake 
spectrum  are  demonstrated  in  the  data  reported  by  Terashima 
(1968)  (see  his  Figs.  5.5  and  6.5).  Unfortunately,  the  defi¬ 
nitions  of  magnitude  and  corner  period  are  different  between 
large  and  small  earthquakes,  and  it  is  impossible  to  decide 
whether  the  apparent  bend  in  the  M  -corner  period  curve  reported 

9 

by  Terashima  is  due  to  the  scale  effect  or  some  other  effects 
such  as  a  gap  in  recording  instrument  response  for  large  and 
small  earthquakes,  or  different  attenuation  and  wave  scattering 
effects  between  short  and  long  periods. 

It  is,  nevertheless,  intriguing  to  consider  the  earthquake 
of  revised  model  A  as  consisting  of  three  distinct  groups: 
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one  with  M  >6,  one  with  6  >  M  >  3,  and  the  other  with 
s  s 

M  <3.  The  largest  earthquake  group  has  roughly  constant 
s 

stress  drop  (10  to  100  bars)  independent  of  earthquake  magni¬ 
tude.  The  medium-sized  group  shares  pre-existing  faults  (a 
few  kilometers  long) ,  and  the  smaller  shock  of  this  group  shows 
less  stress  drop.  In  the  smallest  group,  on  the  other  hand, 
the  stress  drop  increases  with  the  decreasing  magnitude,  in 
agreement  with  Mogi's  (1962)  observation  on  the  size  effect 
on  fracture  strength  and  with  the  high  stress  drop  observed 
in  laboratories  for  small  rock  samples  (the  fracture  strength 
is  proportional  to  L  ^ where  L  is  the  linear  dimension 
of  the  sample,  according  to  Mogi) .  If  the  earth's  crust  con¬ 
tains  weak  zones,  faults  or  cracks  of  a  size  predominantly  a 
few  kilometers,  it  is  conceivable  to  have  such  distinct  groups 
of  earthquakes. 

(2)  Model  B 

Let  us  consider  the  other  extreme  way  of  modifying  the 

w-square  model.  Keeping  the  M  versus  corner-frequency 

s 

curve  unchanged,  and  simply  changing  the  frequency  dependence 
-2  -1 

from  a)  to  w  for  periods  less  than  about  5  seconds,  we 
can  approximately  satisfy  the  observed  Mg  -  relation. 

Since  the  w  ^  dependence  up  to  the  infinite  frequency  results 
in  an  infinite  seismic  energy,  we  assume  that  the  spectral 
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density  drops  sharply  at  about  100  cps.  The  resultant 

family  of  spectral  curves  is  shown  in  Fig.  11. 

This  set  of  curves  can  explain  all  the  main  observations 

discussed  earlier:  (1)  fault  length,  (2)  stress  drop,  (3) 

seismic  moment,  (4)  spectral  ratio  for  T  >10  seconds,  and 

(5)  M  -ft  relation.  In  this  modification,  we  assumed  that 
s  b 

the  suggested  bend  of  M  versus  corner  frequency  may  not  be 

s 

due  to  the  source  effect  but  due  to  instrumental  or  propaga- 

tional  path  effects.  Thus,  the  Mg  -  corner  frequency  relation 

is  the  same  straight  line  as  in  the  w-square  model. 

One  additional  support  of  this  model  comes  from  the  work 

of  Asada  &  Takano  (1963)  and  Takano  (1970,  1971)  on  the 

attenuation  measurement  using  the  spectral  shape  in  the  period 

range  0.1  to  1  seconds.  They  assumed  the  w  1  dependence  at 

the  source  and  obtained  reasonable  Q  values.  They  used  this 

shape  of  source  spectrum  following  Kanai  &  Yoshizawa's  (1958) 

classic  work  on  the  seismic  spectrum  of  nearby  earthquakes 

measured  in  a  deep  mine  in  Japan. 

Instead  of  three  groups  of  earthquakes  for  revised  model 

A,  we  find  two  distinct  groups  for  model  B.  Earthquakes  with 

M  <6  are  all  the  same  kind;  the  w-model  of  Brune  &  King  (1967) 
s 

discussed  earlier.  On  the  other  hand,  earthquakes  with  Mg  >  6 
have  peculiar  spectral  shapes.  For  example,  the  spectrum  for 
M8  first  decreases  as  w  ^  beyond  the  corner  frequency  (about 
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0.003  cps)  but  then  the  decrease  slows  down  to  w  *  beyond 
0.1  cps. 

Within  the  scheme  of  a  simple  dislocation  model  of  an 
earthquake  described  by  a  unidirectional  rupture  propagation 
and  a  step-like  offset  with  a  finite  rise  time  .(Haskell  1964) , 
the  above  w-model  implies  a  very  short  rise  time,  outside  the 
range  of  seismological  observation.  In  view  of  an  incoherent 
rupture  propagation  such  as  described  by  Haskell  (1966)  and 
Aki  (1967) ,  the  co-model  may  correspond  to  a  process  in  which 
the  fault  offset  takes  place  as  a  succession  of  irregular 
rapid  motion  with  a  very  short  time  constant.  In  other  words, 
the  fault  moves  like  a  car  running  at  a  full  speed  on  a  very 
bumpy  road.  Finally,  in  the  framework  of  Brune's  (1970)  model, 
the  w-model  implies  very  small  fractional  stress  drop,  or  very 
large  difference  between  dynamic  and  static  friction.  These 
three  interpretations  appear  to  describe  the  same  phenomena: 
rapid  slips  and  sudden  stops. 

The  peculiar  spectral  shape  for  large  earthquakes  may  be 
explained,  if  the  large  one  behaves  as  an  co-model  up  to  a 
certain  size,  then  transforms  into  an  co-square  model.  It  seems 
reasonable  to  consider  that  rough  and  bumpy  fault  planes  of  the 
co-model  become  smooth  for  larger  displacements  and  produce 
faulting  in  accordance  with  the  co-square  model  which  explains 
the  geologic  observations  on  earthquake  faults. 
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As  mentioned  before#  the  a) -square  model  failed  to  explain 
Tsujiura's  observation  on  two  Aleutian  earthquake  pairs  with 
magnitudes  M  5.7  to  7.0.  The  model  predicted  a  spectral 
ratio  change  of  almost  10  times  from  10  seconds  to  100  seconds, 
while  the  observed  showed  a  very  small  change,  with  a  nearly 
constant  value  over  the  sf-me  period  range .  In  terms  of  the 
revised  model  A,  this  result  means  that  all  these  earthquakes 
probably  had  fault  planes  of  the  same  size.  This  observation 
may  be  approximately  explained  by  the  revised  model  B,  if  the 
transition  from  a  rough  faulting  to  a  smooth  one  takes  place  at 
larger  displacements  or  for  a  larger  fault  size  in  the  Aleutians 
than  in  other  parts  of  the  world. 

Working  models  such  as  A  and  B,  which  predict  the  scaling 
effect  on  seismic  spectrum  over  the  complete  frequency  range 
and  dynamic  range  of  earthquake  seismology,  will  relate  obser¬ 
vations  on  large  earthquakes  with  those  on  small  ones,  and  the 
observations  at  low  frequencies  with  those  at  high  frequencies. 


4.  SCALING  LAW  OF  DISLOCATION  TIME-FUNCTION 

The  source  factor  A(w)  of  the  far -field  seismic  spectrum 
can  be  written  in  terms  of  the  dislocation  D(f,t),  or  the 
displacement  discontinuity,  specified  as  a  function  of  time  t 
and  the  point  f  on  the  fault  plane  £  as  follows  (Haskell  1964) 
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(1) 


where  D(f,u))  is  the  Fourier  transform  of  the  dislocation 
elocity  6(£,t),  r  is  the  distance  to  the  observer  from  the 
surface  element  d^,  and  a  is  the  velocity  of  waves. 

At  low  frequencies,  where  the  change  of  Si£  within  £ 
is  negligible,  the  far-field  spectrum  is  simply  proportional 
to  the  spectrum  of  dislocation  velocity  integrated  over  the 
fault  plane.  Therefore,  the  average  value  of  dislocation 
velocity  spectrum  may  be  obtained  by  dividing  A(u)  by  the 
fault  area. 

For  the  inversion  at  higher  frequencies,  we  shall  assume 
a  simple  model  of  one-dimensional  rupture  propagation,  such  as 
considered  by  Ban-Menahem  (1961) ,  Haskell  (1964)  and  Aki  (1967) . 
In  that  case,  equation  (1)  is  replaced  by 

L  wr 

A(<d)  =  W  D(5,u)  e-i  a  d£ 

J  <2> 

where  K  is  the  coordinate  along  the  path  of  rupture  propa¬ 
gation,  L  is  the  final  fault  length,  W  is  the  width, 

and  dislocation  D(£,w)  j.s  considered  as  the  average  value 
over  the  width. 


For  a  uniform  dislocation  with  constant 


propagation 
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velocity  used  by  Ben-Menahem  (1961) , 


D(£,t)  =  D (t  - 


and  therefore, 


A  ( to )  =  W 


D  (w)  e 


a 


dC 


(3) 


I A  (to)  |=  |  D  (oj)  |  •  WL* 


sin  x 


where  x  = 


La) 


cos  6 


a 


and  6  is  the  angle  between 


the  £-axis  and  the  direction  to  the  observer.  For  high 
frequencies,  therefore,  the  far-field  spectrum  is  propor¬ 
tional  to  the  fault  area,  the  dislocation  velocity  spectrum 
and  a)  .  The  same  result  is  obtained  for  the  incoherent 
propagation  models  used  by  Haskell  and  Aki,  for  frequencies 
beyond  the  corner  frequency.  Thus,  on  the  assumption  of 
one-dimensional  rupture  propagation,  one  can  obtain  the 
spectrum  of  dislocation  velocity  at  frequencies  higher  than 
the  corner  frequency  by  dividing  the  far-field  spectrum  by 
the  fault  area  and  multiplying  it  by  frequency  u). 

The  validity  of  the  above  procedure  should  be  subject 
to  future  critical  investigations.  If  the  rupture  propagation 
is  two-dimensional  as  considered  by  Berckhemer  (1962) , 
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Hirasawa  &  Stauder  (1965)  and  Savage  (1965) ,  one  should 
_2 

expect  w  effect  beyond  the  corner  frequencies.  In  that 

2 

case,  the  correction  should  be  multiplying  w  ,  instead  of 
a).  Of  course,  the  two-dimensional  propagation  is  more 
realistic  because  a  fault  is  not  a  line  but  a  plane.  I  feel 
strongly,  however,  that  one  of  the  corner  frequencies  must 
be  much  higher  than  the  frequencies  associated  with  the 
total  time  of  rupture  propagation  and  the  rise  time  of  fault 
slip,  because  otherwise,  we  expect  that  the  -cube  model  explains 
the  observation  on  the  Mg  -  relation,  but  it  doesn't  (Aki  1967) . 
The  second  corner  frequency  of  rupture  propagation  may  be 
associated  with  a  very  short  transient  time  of  starting 
and  stopping  a  primarily  one-dimensional  rupture  propagation. 

Now  let  us  apply  our  tentative  procedure  of  inversion 
to  the  far-field  spectrum  and  find  the  scaling  law  of  dislo¬ 
cation  time-function.  The  spectral  densities  at  low  fre¬ 
quencies  shown  in  Figs.  1,  9  and  11  are  divided  by  the  fault 
area  (the  square  of  the  corner  period) .  For  frequencies 
higher  than  the  corner  frequency,  an  additional  multiplication 
by  a  factor  proportional  to  w  is  applied.  The  resultant 
family  of  spectral  curves  for  the  dislocation  velocity  spectrum 
are  shown  in  Figs.  12,  13  and  14,  respectively  for  the  u)-squa::e 
model,  revised  model  A  and  B. 

It  is  remarkable  that  the  absolute  value  of  spectral 
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density  of  dislocation  velocity  at  periods  shorter  than  5 

seconds  is  independent  of  magnitude  for  earthquakes  with 

M  >6.5  for  all  three  models, 
s 

According  to  Haskell's  (1969)  calculation  on  the  elastic 
near-field  of  fault  motion,  not  only  the  maximum  displacements 
but  also  the  maximum  velocity  and  acceleration  takes  place  in 
the  immediate  vicinity  of  the  fault  plane.  According  to  Aki's 
(1968)  similar  calculation,  the  seismic  motion  near  a  fault 
depends  neither  on  fault  length  nor  on  fault  width  once  they 
exceed  certain  limits,  but  is  determined  mostly  by  the  dis¬ 
location  time-function  and  velocity  of  rupture  propagation. 

Since  the  rupture  velocity  is  apparently  independent  of 
earthquake  magnitude,  we  must  conclude  that  the  maximum  seismic 
motion  associated  with  an  earthquake  scales  as  the  dislocation 
time-function . 

Then,  our  results  on  the  scaling  law  of  dislocation  time- 
function  will  have  an  important  effect  on  earthquake  engineering, 
because  they  imply  that  the  maximum  seismic  motion  of  an  earth¬ 
quake  in  the  period  range  less  than  5  seconds  is  a  constant 
independent  of  magnitude  for  Mg  >  6.5. 

It  must  be  emphasized  here  that  the  above  conclusion  is 
supported  only  indirectly  by  observations.  No  mention  has 
been  made  of  any  group  of  observations  applying  to  the  upper- 
left  quarter  of  Figs.  1,  9  and  11.  The  waves  in  this  quarter 


I 


-44- 

are  short  waves  coming  from  large  earthquakes.  They  suffer 
not  only  from  the  complexity  of  a  large  source,  but  also 
from  the  complex  path  effect  on  short  waves,  and  this  makes 
the  interpretation  of  the  seismogram  extremely  difficult. 

Because  of  this  difficulty,  they  have  not  attracted  due 
attention.  Housner  (1955)  considered  a  swarm  of  pulses 
random  in  time  resulting  from  release  of  shear  dislocations 
distributed  randomly  over  a  fault  plane.  Haskell  (1964)  intro¬ 
duced  incoherent  rupture  propagation  in  order  to  account  for 
observations  on  short  waves  from  a  large  earthquake.  The  con¬ 
cept  of  a  "multiple  shock”  is  an  old  idea  introduced  to  explain 
this  most  complex  portion  of  a  seismogram  (Stoneley  (1937) , 

Usami  (1956) ,  Wyss  &  Brune  (1967) ,  Trifunac  &  Brune  (1970) , 
among  others) .  Miyamura  et  al.  (1965)  summarize  observations 
and  discuss  physical  mechanisms. 

The  idea  of  a  "multiple  shock"  implies  that  a  large  earth¬ 
quake  consists  of  a  sequence  of  several  small  earthquakes  occur¬ 
ring  within  the  epicentral  region.  One  must  realize,  however, 
the  extremely  large  value  of  seismic  moment  for  the  largest 

earthquakes.  For  example,  the  Alaska  earthquake  of  1964 

29 

(M  «  8.5,  Mq  =  7.5  x  10  dyne  cm)  can  be  a  multiple  of  several 
Rat  Island  earthquakes  of  1965  (M  *  7.9,  MQ  =  1.2  x  102^) ,  but 
requires  about  10,000  San  Fernando  earthquakes  of  1971  (M  =  6.6, 
Mq  =  7.5  x  1025)  .  The  "multiple  shock"  model  of  a  large 
earthquake  will  also  have  a  scale  effect  on  the  spectrum;  the 
spectrum  at  periods  longer  than  the  total  duration  of  faulting 
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will  increase  proportionally  with  the  number  N  of  component 
events,  but  the  high  frequency  spectrum  will  be  proportional 
to  *^N  because  of  random  interference.  Thus,  the  spectral 
ratio  between  large  and  small  earthquakes  at  low  frequency 
will  be  the  square  of  the  ratio  at  high  frequencies.  This 
scale  effect  is  the  same  as  that  of  the  w-model  of  Brune  & 

King  (see  Fig.  5) .  As  discussed  earlier,  observations  for 
large  earthquakes  are  in  favor  of  the  w-square  model,  which 
predicts  that  the  spectral  ratio  at  low  frequency  is  the 
cube  of  the  ratio  at  high  frequency.  The  idea  of  "multiple 
shock"  for  large  earthquakes  is  not  compatible  with  our 
revised  model  B,  in  which  small  earthquakes  are  described  by 
the  w-model,  but  large  ones  by  a  composite  model  of  w  and 
w- square. 
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FIGURE  CAPTIONS 


Source  factor  of  far-field  seismic  spectral  density 
from  earthquakes  with  various  Mg  for  the  w-square 
model  reproduced  from  Aki  (1967) . 

Fault  length  L  as  a  function  of  magnitude  repro¬ 
duced  from  Chinnery  (1969)  with  the  additional 
curve  for  the  w-square  model.  See  text  for  the 
curves  designated  by  the  names  of  investigators. 

The  product  of  fault  length  L  and  the  square  of 
offset  D  as  a  function  of  magnitude  reproduced 
from  King  &  Knopoff  (1968)  ,  with  the  additional 
curve  for  the  w-square  model. 

Seismic  moment  as  a  function  of  magnitude  reproduced 
from  Aki  (1972)  with  additional  lines  (A1)  and  (B) 
for  the  w-model  of  Brune  &  King  (1967)  as  described 

in  text. 

Schematic  representation  of  the  scaling  law  of  seismic 
spectrum  based  on  the  w-model  of  Brune  &  King  (1967)  . 
Locations  of  earthquake  pairs  with  spectral  ratios 
consistent  with  the  w- square  model  (designated  by 
solid  circles)  and  those  inconsistent  (designated  by 


crosses) . 
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Comparison  of  observed  Mg  -  relation  with 

the  one  predicted  for  the  w-square  model. 

8  Theoretical  Mg  -  ^  relation  for  the  w-square 
model  as  compared  with  empirical  formulas  of 
Gutenberg-Richter  (1956)  and  Evernden  et  al.  (1971) 
for  earthquakes  and  that  of  Thirlaway-Carpenter 
(1966)  for  explosions. 

9  Source  factor  of  far-field  seismic  spectral  density 
from  earthquakes  with  various  Mg  for  the  revised 
model  A. 

10  Theoretical  M  versus  corner- frequency  relation 

s 

for  the  w-square  model  as  compared  with  those 
implied  by  the  Mg  -  relations  of  Gutenberg- 
Richter  (1956)  and  Evernden  et  al.  (1971) . 

11  Source  factor  of  far-field  seismic  spectral  density 
from  earthquakes  with  various  Mg  for  the  revised 
model  B. 

12  Spectral  densities  of  dislocation  velocity  for 

different  M  based  on  the  w-square  model, 
s 

13  Spectral  densities  of  dislocation  velocity  for 

different  M  based  on  the  revised  model  A. 
s 

Spectral  densities  of  dislocation  velocity  for 

different  M  based  on  the  revised  model  B. 
s 


FIGURE  14 


REFERENCES 


Abe,  K.,  1970.  Determination  of  seismic  moment  and  energy 
from  the  earth's  free  oscillation,  Phys.  Earth  Planet. 
Interiors,  4,  49-61. 

Aki,  K.,  1966.  Generation  and  propagation  of  G  waves  from  the 
Niigata  earthquake  of  June  16,  1964,  Bull.  Earthq.  Res. 
Inst. ,  Univ.  Tokyo,  £4,  23-88. 

Aki,  K.,  1967.  Scaling  law  of  seismic  spectrum,  J.  Geophys. 
Res. ,  72,  1217-1231. 

Aki,  K.,  1968.  Seismic  displacements  near  a  fault,  J .  Geophys . 
Res.,  73,  5359-5376. 

Aki,  K.,  1969.  Analysis  of  the  seismic  coda  of  local  earth¬ 
quakes  as  scattered  waves,  J .  Geophys .  Res . ,  74,  615-631. 
Aki,  K.,  1972.  Earthquake  Mechanism,  Proceeding  of  the  Final 
UMP  Symposium,  Moscow  1971,  ed.  by  R.  Ritsema. 

Archambeau,  C.,  1968.  General  theory  of  elastodynamic  source 
fields.  Rev.  Geophys.,  £,  241-288. 

Asada,  T.,  1957.  Observations  of  nearby  microearthquakes  with 
ultra-sensitive  seismometers,  J.  Phys.  Earth,  5,  03-113. 
Asada,  T.  &  Takano,  K.,  1963.  Attenuation  of  short-period  P 
waves  in  the  mantle,  J.  Phys.  Earth,  11 ,  25-34. 

Basham,  P.W. ,  1969.  Canadian  magnitudes  of  earthquakes  and 
nuclear  explosions  in  south-western  North  America, 


-50- 

Bath,  M.,  &  Duda,  S.J.,  1964.  Earthquake  volume,  fault  plane 
area,  seismic  energy,  strain,  deformation  and  related 
quantities,  Ann.  Geofis.  Rome,  17 ,  353-368. 

Benioff,  H.,  1951.  Earthquakes  and  rock  creep.  Part  I:  creep 
characteristics  of  rocks  and  the  origin  of  aftershocks, 

Bull.  Seismol.  Soc.  Am.,  41,  31-62. 

Ben-Menahem,  A.,  1961.  Radiation  of  seismic  surface  waves  from 
finite  moving  sources,  Bull.  Seismol.  Soc.  Am.,  51,  401-435. 

Ben-Menahem,  A.,  &  Toksfiz,  M.N.,  1963.  Source  mechanism  from 
spectrum  of  long-period  surface  waves:  2.  Kamachatka 
earthquake  of  November  5,  1952,  J.  Geophys.  Res.,  6jJ, 
5207-5222. 

Berckhemer,  H.,  1962.  Die  Ausdehnung  der  BruchflSche  im 

Erdbebenherd  und  ihr  Einfluss  auf  das  seismiche  Wellen- 
spektrum,  Gerlands  Beitr.  Geophys.,  71,  5-26. 

Brune,  J.N.,  1968.  Seismic  moment,  seismicity,  and  rate  of 
slip  along  major  fault  zones,  J.  Geophys.  Res., 

777-784. 

Brune,  J.N.,  1970.  Tectonic  stress  and  the  spectra  of  seismic 
shear  waves  from  earthquakes,  J.  Geophys.  Res.,  7j>, 
4997-5009. 

Brune,  J.N.,  1971.  Seismic  sources,  fault  plane  studies  and 
tectonics,  Trans.  Am.  Geophys.  Union,  52 ,  178-187. 

Brune,  J.N.,  Correction  to  "Tectonic  stress  and  the  spectra 

of  seismic  shear  waves  from  earthquakes,"  J.  Geophys.  Res., 
76,  5002. 


£aseat*(fai 


t 


-51- 

Brune,  J.N.,  &  Allen,  C.R.,  1967.  A  low-stress-drop,  low- 

magnitude  earthquake  with  surface  faulting:  The  Imperial, 
California,  earthquake  of  March  4,  1966,  Bull.  Seismol. 

Soc.  Am.,  57 ,  501-514. 

Brune,  J.N.,  &  Engen,  G.,  1969.  Excitation  of  mantle  Love 

waves  and  definition  of  mantle  wave  magnitude,  Bull.  Seismol. 
Soc.  Am.,  59,  923-933. 

Brune,  J.N.,  &  King,  D.Y.,  1967.  Excitation  of  mantle  Rayleigh 
waves  of  period  100  seconds  as  a  function  of  magnitude, 

Bull.  Seismol.  Soc.  Am.,  57 ,  1355-1365. 

Capon,  J. ,  Greenfield,  R.J.  &  Lacoss,  R.T.,  1967.  Surface- 
versus  body-wave  magnitude  results,  in  Semiannual 
Technical  Summary,  Lincoln  Laboratories,  M.I.T.,  pp.  3-5, 
Lexington,  Mass.,  June  30. 

Chalturin,  V.I.,  Sept.  1970.  Proceeding  of  Eurpopean  Eng. 
Seismological  Commission,  Sofia,  Bulgaria. 

Chinnery,  M. ,  1964.  The  strength  of  the  earth's  crust  under 

horizontal  shear  stress,  J.  Geophys.  Res.,  6f),  2085-2089. 

Douglas,  B.M. ,  &  Ryall,  A.,  1972.  Spectral  characteristics 

and  stress  drop  for  microearthquakes  near  Fairview  Peak, 
Nevada,  J.  Geophys.  Res.,  77 ,  351-359. 

Evernden,  J.F.,  Best,  W.J.,  Pomeroy,  P.W.,  McEvilly,  T.V., 

Savino,  J.M.,  &  Sykes,  L.R. ,  1971.  Discrimination  between 
small-magnitude  earthquakes  and  explosions,  J.  Geophys.  Res., 
76,  8042-8055. 


-52- 

Filson,  J.,  &  McEvilly ,  T.V. ,  1967.  Love  wave  spectra  and  the 
mechanism  of  the  1966  Parkfield  sequence,  Bull.  Seismol. 

Soc  .  Am. ,  !57 ,  1245-1257  . 

Furuya,  I., 1969.  Predominant  period  and  magnitude,  J.  Phys. 

Earth,  17,  119-126. 

Gutenberg,  B.,  &  Richter,  C.F.,  1956.  Earthquake  magnitude, 

intensity,  energy,  and  acceleration,  2,  Bull.  Seismol.  Soc. 

Am. ,  46,  105-145. 

Haskell,  N.,  1964.  Total  energy  and  energy  spectral  density 
of  elastic  wave  radiation  from  propagating  faults,  Bull . 

Seismol.  Soc.  Am.,  56,  1811-1842. 

Haskell,  N.,  1966,  Total  energy  and  energy  spectral  density  of 

elastic  wave  radiation  from  propagating  faults,  2,  A  statis¬ 
tical  source  model,  Bull.  Seismol.  Soc.  Am.,  56 ,  125-140. 

Haskell,  N.,  1969.  Elastic  displacements  in  the  near-field  of 
a  propagating  fault.  Bull.  Seismol.  Soc.  Am.,  59 ,  865-908. 
Hirasawa,  T.,  &  Stauder,  W. ,  1965.  On  the  seismic  body  waves 
from  a  finite  moving  source ,  Bull.  Seismol.  Soc.  Am. ,  55 , 

237-262. 

Housner,  G.W.,  1955.  Properties  of  strong  ground  motion  earth¬ 
quakes,  Bull.  Seismol.  Soc.  Am.,  45 ,  197-218. 

Ida,  Y.,  &  Aki,  K.,  1972.  Seismic  source  time  function  of  propa¬ 
gating  longitudinal-shear  cracks,  J.  Geophys.  Res.,  76  (in  press). 
Iida,  K.,  1959.  Earthquake  energy  and  earthquake  fault,  J. 

Earth  Sci.,  Nagoya  Univ.,  7,  98-107. 

Iida,  K.,  1965.  Earthquake  magnitude,  earthquake  fault,  and 

source  dimensions,  J.  Earth  Sci.,  Nagoya  Univ.,  13^,  115-132. 


1 


-53- 


Kanai,  K. ,  &  Yoshizawa,  S.,  1958.  The  amplitude  and  the  period 

of  earthquake  motions,  Bull.  Earthq.  Res.  Inst.,  3£,  275-293. 

Kanamcri,  H.,  1970a.  Synthesis  of  long-period  surface  waves  and 
its  application  to  earthquake  source  studies  -  Kurile 
Islands  earthquake  of  October  13,  1963.  J.  Geophys.  Res., 

75,  5011-5027. 

Kanamori,  H.,  1970b.  The  Alaska  earthquake  of  1964:  Radiation 
of  long-period  surface  waves  and  source  mechanism,  J. 

Geophys.  Res.,  75 ,  5029-5040. 

Kanamori,  H.,  1971.  Seismological  evidence  for  a  lithospheric 
normal  faulting  -  the  Sanriku  earthquake  of  1933,  Ph^s. 

Earth  Planet.  Interiors,  4,  289-300. 

King,  C.Y.,  &  Knopoff,  L.,  1968.  Stress  drop  in  earthquakes, 

Bull.  Seismol.  Soc.  Am.,  58 ,  249-257. 

Lieberman,  R.C.,  &  Basham,  P.W.,  1971.  Excitation  of  surface 

waves  by  the  Aleutian  underground  explosion  Milrow  (October 
2,  1969),  J.  Geophys.  Res.,  76 ,  4030-4034. 

Lieberman,  R.C. ,  King,  C.Y.,  Brune,  J.N.,  &  Pomeroy,  P.W. ,  1966. 
Excitation  of  surface  waves  by  the  underground  nuclear 
explosion  Long  Shot,  J.  Geophys.  Res.,  7_1,  4333-4339. 

Lieberman,  R.C.,  &  Pomeroy,  P.W.,  1967.  Excitation  of  surface 

waves  by  events  in  southern  Algeria,  Science,  156 ,  1098-1100. 

Lieberman,  R.C.,  &  Pomeroy,  P.W.,  1969.  Relative  excitation  of 
surface  waves  by  earthquakes  and  underground  explosions, 

J.  Geophys.  Res.,  74 ,  1575-1590. 

Lieberman,  R.C.,  &  Pomeroy,  P.W.,  3970.  Source  dimensions  of 


/ 


-54- 


small  earthquakes  as  determined  from  the  size  of  after¬ 
shock  zone,  Bull.  Seismol.  Soc.  Am.,  60 ,  879-896. 

Linde,  A.T.,  &  Sacks,  I.S.,  1971.  Errors  in  the  spectral 

analysis  of  long-period  seismic  body  waves,  J.  Geophys. 

Res. ,  76,  3226-3336. 

Miyamura,  S.,  Omote,  S.,  Teisseyre,  R. ,  &  Vesanen,  E.,  1965. 
Multiple  shocks  and  earthquake  series  pattern,  Bull. 

Intern.  Inst.  Seis.  Earthq.  Eng.,  2,  71-92. 

Mogi,  K. ,  1962.  The  influence  of  the  dimensions  of  specimens 
on  the  fracture  strength  of  rocks  -  comparison  between 
the  strength  of  rock  specimens  and  that  of  the  earth's 
crust,  Bull.  Earthq.  Res.  Inst.,  Univ.  Tokyo,  40^,  175-185. 

Molnar,  P.,  Savino,  J.,  Sykes,  L.R. ,  Lieberman,  R.C.,  Hade,  G. , 

&  Pomeroy,  P.W.,  19  .  Small  earthquakes  and  explosions 

in  western  North  America  recorded  by  new  high-gain,  long- 
period  seismographs,  Nature ,  224 ,  1268-1273. 

Otsuka,  M. ,  1965.  Earthquake  magnitude  and  surface  fault 
formation,  J.  Seis.  Soc.  Japan,  Ser.  2,  18,  1-8. 

Press,  F.,  1967.  Dimensions  of  the  source  region  for  small 

shallow  earthquakes,  Proceedings  of  the  VESIAC  Conference 
on  the  Current  Status  and  Future  Progress  for  Understanding 
the  Source  Mechanism  of  Shallow  Seismic  Events  in  the  3 
to  5  Magnitude  Range. 

Savage,  J.C.,  1966.  Radiation  from  a  realistic  model  of 
faulting,  Bull.  Seismol.  Soc.  Am.,  56,  577-592. 


-55- 


Stauder ,  W. ,  1968.  Tensional  character  of  earthquake  foci 
beneath  the  Aleutian  trench  with  relation  to  sea-floor 
spreading,  J.  Geophys.  Res..  73i,  7693-7701. 

Stoneley,  R. ,  1J37.  The  Mongolian  earthquake  of  1931,  August 
10.  Brit.  Ass.  Seis.  Conun.  42nd  Rep.,  5-6. 

Takano,  K.,  1970.  "Attenuation  of  short  period  seismic  wave3 
in  the  upper  mantle  and  its  regional  difference,  J .  Phys . 
Earth,  18,  171-180. 

Takano,  K. ,  1971.  A  note  on  the  attenuation  of  short  period  P 
and  S  waves  in  the  mantle,  J.  Phys.  Earth.  19,  155-164. 

Takano,  K. ,  1971.  Analysis  of  seismic  coda  waves  of  ultra 
microearthquakes  in  the  Matsushiro  area  -  a  comparison 
with  Parkfield ,  California,  J.  Phys.  Earth,  19,  209-215. 

Terashima,  T.,  1968.  Magnitude  of  microearthquake  and  the 

spectra  of  microearthquake  waves,  Bull.  Intern.  Inst.  Seism. 
Earthq.  Eng.,  5,  31-108. 

Thirlaway,  H.I.S.,  &  Carpenter,  E.W. ,  1966.  Seismic  signal 
anomalies,  travel  times,  amplitudes,  and  pulse  shapes. 
Proceedings  of  the  VESIAC  Special  Study  Conference  on 
Seismic  Signal  Anomalies,  Travel  Times,  Amplitudes,  and 
Pulse  Shapes,  Beaugency,  France,  October. 

Tocher,  D.,  1958.  Earthquake  energy  and  ground  breakage.  Bull. 
Seismol .  Soc.  Am..  48,  147-152. 


-56- 


Tocher,  D.,  1960.  Movement  on  faults,  Proc.  2nd  World  Conf. 
Earthquake  Engineering,  1,  551-564. 

Trifunac,  M. ,  &  Brune,  J.N.,  1970.  Complexity  of  energy  re¬ 
lease  during  the  Imperial  Valley,  California,  earthquake 
of  1940,  Bull.  Seismol.  Soc.  Am.,  (SO,  137-160. 

Tsai,  Y.B.,  &  Aki,  K. ,  1969.  Simultaneous  determination  of 
the  seismic  moment  and  attenuation  of  seismic  surface 
waves.  Bull.  Seismol .  Soc.  Am.,  59,  275-287. 

Tsai,  Y.B.,  &  Aki,  K.,  1970a.  Source  mechanism  of  the  Truckee, 
California  earthquake  of  September  12,  1966,  Bull.  Seismol. 
Soc.  Am. ,  60,  1199-1208. 

Tsai,  Y.B.,  &  Aki,  K.,  1970b.  Precise  focal  depth  determination 
from  amplitude  spectra  of  surface  waves,  J .  Geophys .  Res . , 
75,  5729-5743. 

Tsai,  Y.B.,  &  Aki,  K.,  1971.  Amplitude  spectra  of  surface  waves 
from  small  earthquakes  and  underground  nuclear  explosions, 
J.  Geophys.  Res.,  76 ,  39  0-3952. 

Tsuboi,  C.,  1940.  Isostasy  and  maximum  earthquake  energy,  Proc. 
Imp.  Acad . ,  16,  449. 

Tsuboi,  C.,  1956.  Earthquake  energy,  earthquake  volume,  after¬ 
shock  area  and  strength  of  the  earth's  crust,  J .  Phys . 
Earth,  £,  63. 

Tsuboi,  C.,  1958.  On  seismic  activities  in  and  near  Japan, 

"Contributions  in  Geophysics;  in  honor  of  Beno  Gutenberg, " 
Pergamon  Press. 


-57- 


Tsuboi,  C.,  1965.  A  theoretical  derivation  of  the  magnitude- 
energy  equation  of  earthquakes,  log  E  =  aM  +  b,  Proc. 

Japan  Acad . ,  41,  588. 

s 

Tsujiura,  M. ,  1966.  Frequency  analysis  of  seismic  waves  (1), 

Bull.  Earthq.  Res.  Inst.,  Tokyo  Univ.,  4£,  873-891. 

Tsujiura,  M. ,  1967.  Frequency  analysis  of  seismic  waves  (2), 
Bull.  Earthq.  Res.  Inst..  Tokyo  Univ.,  45,  973-995. 

Tsujiura,  M. ,  1969.  Regional  variation  of  P  wave  spectra  (1), 
Bull.  Earthq.  Res.  Inst.,  Tokyo  Univ.,  47,  613-633. 

Tsujiura,  M. ,  1972.  Spectrum  of  seismic  waves  and  its  dependence 
on  magnitude  (1),  Bull.  Earthq.  Res.  Inst.,  Tokyo  Univ. 

(in  press) . 

Usami,  T. ,  1956.  Seismometrical  study  of  Boso-oki  earthquake 
of  November  26,  1953,  Quart.  J.  Seismology  (Kenshin-Jiho) 
Japan  Meteor.  Agency,  21,  3. 

Wyss,  M.,  1970.  Observation  and  interpretation  of  tectonic 
strain  release  mechanisms,  Ph.D.  thesis,  Ca-if.  Inst, 
of  Technology. 

Wyss,  M. ,  &  Brune,  J.N.,  1967.  The  Alaska  earthquake  of  28 
March  1964:  A  complex  multiple  rupture,  Bull.  Seismol. 

Soc.  Am.,  57,  1017-1023. 

Wyss,  M. ,  &  Brune,  J.N.,  1968.  Seismic  moment,  stress,  and 
source  dimensions  for  earthquakes  in  the  California- 
Nevada  region,  J.  Geophys.  Res.,  73 ,  4681-4694. 


SEISMIC  MOMENT  IN  DYNE-CM 


a  p  n  4. !  ubD|/\| 


log  (LD  )(  L,  D  in  cm) 


SEISMIC  MOMENT  IN  DY N E  CM 


Fig.  4 


1.  Central  America 

2.  North  Atlantic  7  / 

3.  Mexico -Guatemala  /  /  \  // 

4  Peru  \  H 

(Evernden  et  aU  /  f  S/ /  \ \  \  l 

/i//Ay  w 

yjk  A]  / 

//m  /  /  / 

//  //  //  \  /v 

;  //by  kt 

/  7  fit 

1  /  f  i/  /  #  / 

/  V  if  /  /I  / 

Model 


4  Peru 


5.  U.  S.  (Basham,  Everden  et  al,  Lieberman- 
Pomeroy) 


6.  Aleutian-  Kamchatta 
7  Central  Asia 


(Lieberman-  Pomeroy) 


Ms  defined 


REVISED 
MODEL  "A" 


?IOD  IN  SEC 


(  33V0S  3AllV13a)  A1ID013A 
NOIIVOOISIQ  JO  A1ISN3G  !Vai03dS 


Fig.  12 


•mu*/**’*"' 


PERIOD  IN  SEC 


-70- 


o  o  - 

o  - 

(31V0S  3  Ai±V13d  )  A1ID013A 
N0I1VD01SIQ  dO  A1ISN3G  “IVdlD3dS 


Fig.  13 


PERIOD  IN  SEC 


PERIOD  IN  SEC 


.  jpmnpt — -  - - - 

-72- 

3 . 6  Tectonic  Strain  Release  by  Underground  Nuclear  Explosions 
and  its  Effect  on  Seismic  Discrimination  by  M.  Nafi 
Toksttz  and  Harold  H.  Kehrer 

Summary 

The  analysis  of  surface  waves  from  a  large  number  of 
underground  nuclear  explosions  reinforces  the  hypothesis  of 
tectonic  strain  release.  Such  strain  release  is  dependent 
on  rock  type  and  ambient  stress  levels.  In  harder  media, 
such  as  granite,  strain  release  is  considerably  greater  than 
in  soft  media  such  as  loose  alluvium  and  salt.  In  the  case 
of  three  events  for  which  the  tectonic  component  exceeded 
that  of  the  explosion  itself,  the  effect  on  the  Mg  -  m^ 
discriminant  due  to  the  added  surface  wave  energy  was  not 
noticeable.  This  was  probably  due  to  the  averaging  of 
magnitude  over  all  azimuths.  A  decaying  pulse-type  source¬ 
time  function  is  consistent  with  surface  wave  amplitude 
spectra  and  can  explain  the  long-period/short-period  spectral 
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INTRODUCTION 

The  source  mechanisms  of  underground  nuclear  explosions 
have  been  of  great  seismological  interest  since  it  was  first 
observed  that  a  pure  explosive  source  alone  could  not  adequately 
account  for  the  radiation  patterns  of  seismic  waves  recorded 
at  distant  stations.  Most  notable  among  these  observations 
are  the  presence  of  SH  and  Love  waves  and  azimuthal  asymmetry 
of  the  Rayleigh  wave  radiation  patterns.  There  have  been  a  num¬ 
ber  of  studies  dealing  with  the  radiation  patterns  and  the 
generation  of  SH-type  seismic  waves  by  underground  nuclear 
explosions  (Press  &  Archambeau  1962;  Brune  &  Pomeroy 
1963;  Aki  1964;  Toksdz  et  al,  1964,  1965;  Toksoz  1967;  Archam¬ 
beau  1968;  Kehrer  1969;  Molnar  et  al  1969;  Archambeau  & 

Sammis  1970;  Toksdz  &  Kehrer  1971;  Tcksdz  et  al  1971)  .  It  is 
generally  accepted  that  the  Love  waves  are  generated  at  or 
very  near  the  source  and  that  they  may  be  associated  with  the 
release  of  tectonic  strain  energy  by  the  explosion.  As  recent 
calculations  by  Leavy  (1971)  show,  near-source  inhomogeneities 
and  scattering  cannot  account  for  the  observed  amplitudes  of 
Love  waves  Geological  studies  (McKeown  et  al  1966;  Barosch 
1968;  Bucknam  1969;  McKeown  &  Dickey  1969;  Dickey  1968,  1969, 
1971)  clearly  indicate  fault  movements  and  sur¬ 
face  displacements  associated  with  the  majority  of  explosions. 
Aftershocks  that  follow  the  explosions  (Hamilton  &  Healy 
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1969;  Ryall  &  Savage  1969;  Stauder  1971)  generally  resemble 
the  sequence  that  follows  earthquakes  and  they  indicate  the 
disturbance  of  the  regional  tectonic  stress  patterns  by  the 
explosion. 

When  observed  outside  the  source  region,  seismic  waves 
which  are  due  to  a  tectonic  strain  release  component  will  have 
spectral  properties  and  radiation  patterns  similar  to  those 
of  an  earthquake.  These  will  be  superimposed  on  waves  radiated 
isotropically  by  the  pure  explosive  source.  The  questions 
that  arise  are:  1)  what  is  the  extent  of  tectonic  strain 
energy  release  and  2)  how  significantly  does  this  complication 
affect  some  of  the  seismic  discriminants  such  as  Ms  -  m^ 
relationships?  In  this  paper  we  address  ourselves  to  these  points. 

! 

I 

The  paper  is  divided  into  three  parts.  In  Part  I,  we 
present  the  strain  release  characteristics  of  about  20  under¬ 
ground  explosions.  For  this,  we  use  Rayleigh  wave  radiation 
patterns  as  well  as  Love-to-Rayleigh  wave  amplitude  ratios.  We 
show  the  importance  of  medium  properties  and  tectonic  framework 
in  controlling  the  energy  release.  In  Part  II,  we  evaluate 
the  effect  of  the  released  strain  energy  on  the  surface  wave 
magnitudes.  Part  III  is  devoted  to  a  somewhat  separate  subject  — 
the  source  function  of  the  explosion.  A  decaying  pulse-type 
pressure  function  fits  the  observed  surface  spectra  and  explains 
the  success  of  some  of  the  spectral  ratio  discriminants. 

I 
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I.  TECTONIC  STRAIN  RELEASE  CHARACTERISTICS  OF  EXPLOSIONS 

Examples  of  Rayleigh  and  Love  waves  generated  by  three 
explosions  are  shown  in  Fig.  1.  These  illustrate  the  differ¬ 
ences  in  the  relative  generation  of  Love  waves  by  NTS  (Nevada 
Test  Site)  explosions  recorded  at  Weston,  Massachusetts.  Love 
waves  appear  on  the  tangential  components  (T) ,  about  two  minutes 
earlier  than  the  Rayleigh  waves  on  the  vertical  components  (2) . 
The  large  amplitude  of  Love  waves  relative  to  Rayleigh  waves 
for  Pile  Driver  and  Greeley  indicates  the  greater  amount  of 
tectonic  strain  release  by  these  explosions. 

To  determine  the  contribution  of  strain  release  on  the 
radiated  wave,  we  formulate  a  composite  source  consisting  of 
a  pure  explosion  plus  a  double-couple  component  representing 
the  strain  release.  For  such  a  source,  the  theoretical  formu¬ 
lations  for  the  azimuthal  dependence  of  Rayleigh  wave  amplitudes 
and  Love/Rayleigh  amplitude  ratios  can  be  written  as  (Toksflz 
et  al  1965,  1971;  Kehrer  1969): 
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where  A  and  Ar  are  the  Rayleigh  and  Love  amplitude  factors 
R  l 

(Harkrider  1964);  kR  and  kL  are  the  Rayleigh  and  Love 

wave  numbers;  and  are  components  of  particle  velocity 

at  the  surface;  aQ  is  Poisson's  ratio  at  the  surface; 
and  yT  are  the  Rayleigh  and  Love  wave  attenuation  coefficients; 
6  and  X  are  the  dip  and  slip  directions  of  the  fault  plane; 

9  is  its  azimuth;  and  F  is  the  strength  of  the  double-couple 
relative  to  the  explosion.  In  Eq.  (1) »  the  difference  in  ex¬ 
plosion  and  earthquake  source  time  functions  and  source  dimen¬ 
sions  are  considered  negligible. 

For  a  horizontal  double-couple,  6  =  90°  and  X  -  0°  and 
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In  this  study  we  use  Eq.  (2)  to  determine  the  strength  of  the 
double-couple  (F  value)  and  the  orientation  of  the  double¬ 
couple  (  0  )  from  Rayleigh  wave  amplitudes  and  Love/Ray  i.eigh 
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amplitude  ratios.  Other  quantities  are  computed  theoretically 
for  specified  layered  earth  models  (Harkrider,  1964) 

Long-period  seismograms  were  obtained  from  stations  through¬ 
out  North  America  for  the  study  of  the  U.  S.  explosions.  These 
included  stations  of  the  World  Wide  Standard  Seismograph  Network  (WWSSN)  , 
the  Long  Range  Seismic  Measurements  Network  (LRSM) ,  the  Canadian 
Standard  seismograph  Network  (CSSN) ,  and  two  additional  stations  at 
Berkeley  and  Pasadena  which  provided  long-period  data.  Figure 
2  shows  the  distribution  of  recording  stations  in  North  America. 
However ,  not  all  stations  provided  useful  data  for  each  event. 

For  the  presumed  Soviet  events,  data  was  obtained  from  the  WWSSN 
stations  in  Europe  and  Asia. 

The  film  records  of  those  events,  for  which  the  surface 
wave  data  were  of  sufficient  quality,  were  digitized,  band-pass 
filtered  (8-60  sec)  to  reduce  noise,  and  the  two  horizontal  com¬ 
ponents  were  rotated  to  radial  and  tangential  directions  with 
respect  to  the  epicenter.  In  this  way  essentially  pure  Love 
waves  were  obtained  on  the  tangential  component.  Figure  3  is 
an  example  of  this  analysis  for  the  explosion  Boxcar  recorded 
at  Mould  Bay,  Canada.  The  three  lower  traces  are  plotted  from 
digital  data  and  show  the  unfiltered  north-south,  east-west, 
and  vertical  components.  The  three  upper  traces  show  the  fil¬ 
tered  data,  where  the  horizontal  components  have  been  rotated 
to  radial  and  tangential  directions. 
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The  tangential  and  the  vertical  components  were  then 
Fourier-analyzed  to  yield  the  amplitude  spectra.  The  proper 
time  windows  for  this  analysis  were  determined  by  examining 
the  filtered  trace  at  each  station  individually  and  noting  the 
onset  and  duration  of  the  wave  trains. 

Love  wave  to  Rayleigh  wave  amplitude  ratios  were  deter¬ 
mined  by  computing  the  spectral  ratio  of  the  tangential  com¬ 
ponent  to  the  vertical  component  over  the  period  range  8  to 
60  seconds.  In  general,  the  spectral  ratios  were  somewhat 
oscillatory  and  thus  the  ratio  at  any  one  particular  period 
could  not  be  expected  to  be  reliable  at  all  stations.  The 
period  range  was,  therefore,  divided  into  three  bands,  from 
9  to  15  sec,  from  15  to  22  sec,  and  from  22  to 
30  sec.  The  average  value  of  the  ratio  in  each  band  was  then 
taken  to  represent  the  ratio  over  that  band. 

An  automatic  error  scheme  (Kehrer  1969)  was  used  to  deter¬ 
mine  the  best  fit  of  Eq.  (2)  to  the  observed  data.  In  applying 

Eq.  (2),  the  term  k^2AL/k^2AR (~)  was  taken  to  be  a  constant 

o 

at  all  stations.  An  error  term  E  was  then  formed  between  the 
observed  and  the  theoretical  ratios  for  combinations  of  F  and 
6  . 


n 

l 

i=l 


S'Fj  cos  20^ 

1  +  Fj  sin  20^ 


./  N 


E. 

3 


1/2 


(3) 


-80- 


The  jth  combination  of  F  and  0  which  fits  the  data  best 
will  be  that  which  minimizes  E.  Lj/Ri  t^le  experimentally 
measured  Love  wave  to  Rayleigh  wave  amplitude  ratio  at  a  par¬ 
ticular  station  i,  S  is  a  constant,  and  N  is  the  number 
of  stations.  The  values  of  E  were  contoured  on  a  Stromberg- 
Carlson  4020  grid  for  the  three  period-bands  of  each  event. 

Figure  4  shows  these  results  for  the  explosion  Faultless.  The 
minima  indicate  the  orientation  of  the  best  fitting  right-lateral, 
vertical,  strike-slip  fault  as  measured  clockwise  from  the  north 
at  the  source.  A  comparison  of  the  error  plots  over  the  three 
bands  in  Fig.  4  reveals  the  frequency  dependence  of  F.  The 
orientation  0  6f  the  double-couple  remains  essentially  un¬ 
changed  over  the  three  bands  while  F  is  greatest  at  long 
periods  and  smallest  at  short  periods.  This  relationship  was 
observed  for  nearly  all  events  and  indicates  that  the  source 
of  Love  waves  is  less  efficient  at  short  periods  than  the  Ray¬ 
leigh  wave  source.  A  similar  phenomenon  was  noted  by  Aki  et 
al  (1969)  for  body  waves  from  the  Benham  event.  Figure  5  shows 
the  relationship  oetween  F  and  period  for  several  events. 

Observed  and  calculated  Love/Rayleigh  amplitude  ratios 
as  a  function  of  azimuth  for  three  explosions  are  shown  in  Fig. 

6.  A  summary  of  the  best  fitting  F  and  0  in  the  15  to  22 
sec  period  range  for  each  event  is  given  in  Table  1.  For  some 
explosions  (Buff  and  the  presumed  Soviet  events  of  3-20-66, 
3-26-67,  and  10-21-67)  for  which  the  surface  wave  data  was  not 
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of  sufficient  quality  to  permit  the  calculation  of  spectral 
ratios,  the  parameters  F  and  0  were  determined  from  mea¬ 
surements  of  the  maximum  amplitudes  of  Love  and  Rayleigh  waves. 
Included  in  Table  1  is  the  ratio  of  the  surface  wave  energy 
of  the  double-couple  to  the  energy  of  the  explosion.  This 
ratio  (Toksfiz  etal  1965)  is  closely  approximated  by 

E.  ./E0  =  4/3  Fz.  (4) 

tect  exp 

i 

Several  conclusions  can  be  drawn  from  Table  1.  The  close 

i 

agreement  in  fault  orientation  of  most  of  the  events  at  the 
Nevada  Test  Site  indicates  that  a  regional  stress  field  is 
probably  controlling  the  strain  release.  The  large  variation 
within  the  two  Soviet  testing  areas  may  be  due  to  a  lack  of  sufficien 

I 

data  in  the  case  of  three  of  the  events.  The  dependence  of  F 
and  hence  the  amount  of  stress  release,  on  lithology  was  pre¬ 
viously  noted  by  Toksflz  (1967)  and  by  Toksdz  &  Kehrer  (1971). 

Table  1  provides  additional  evidence  that  the  multipolar  source 
component  is  larger  in  more  competent  rock.  Thus  strain  re¬ 
lease  from  a  particular  explosion  is  dependent  on  the  rock 
medium  and  the  regional  state  of  strain.  It  is  also  apparent 
that  there  are  variations  for  a  given  rock  type,  as  in  the  case 
of  the  three  events  in  granite.  This  can  be  explained  by  the 
differences  in  ambient  stress  levels  and  available  strain  energy 
(Toksdz  et  al  1971).  When  the  F  value  exceeds  unity  the 
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Rayleigh  wave  radiation  pattern  will  show  a  lobate  structure 
and  reversed  polarity  in  alternating  quadrants.  This  has  been 
discussed  in  an  earlier  paper  (Toksdz  &  Kehrer  1971)  and  is 
here  demonstrated  in  Fig.  7  for  the  Pile  Driver  and  Greeley 
explosions. 

The  north-northwest  orientation  of  the  fault  planes  of 
most  of  the  Nevada  Test  Site  explosions  is  in  good  agreement 
with  lineaments  and  dominant  fault  patterns  in  the  area.  Figure 
8  shows  the  Yucca  Flat  portion  of  the  Nevada  Test  Site  with  the 
source  mechanisms  of  the  explosions  studied  there.  Included 
are  the  major  natural  faults  and  some  of  the  fractures  produced 
by  the  events.  The  chief  feature  of  Yucca  Flat  is  the  Yucca 
Fault  which  trends  northward  through  alluvium  and  is  of  recent 
age.  Although  most  observed  displacements  across  the  Yucca 
Fault  are  vertical  following  nearby  explosions,  numerous  en 
echelon  fractures  along  the  fault  were  produced  in  the  alluvium 
by  Corduroy  (Barosch  1968).  Such  an  en  echelon  pattern  is  gen¬ 
erally  related  to  strike-slip  movement  in  the  underlying  base¬ 
ment.  The  trend  of  the  pattern  in  this  case  suggests  right- 
lateral  slippage. 

Figure  9  shows  the  Pahute  Mesa  portion  of  the  Nevada  Test 
Site  with  the  natural  faults,  explosion-produced  fractures,  and 
the  events  studied.  The  main  structural  feature  of  Pahute  Mesa 
is  the  Silent  Canyon  Caldera,  which  encloses  the  five  events. 
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Many  normal  faults,  striking  north-northwest,  cut  the  thick 
sequence  of  Tertiary  volcanic  rock.  Recent  natural  movement 
along  some  of  these  faults  has  been  inferred  (McKeown  et  al 
1966).  For  Benham  and  Boxcar,  studies  of  aftershock  distribu¬ 
tion  (Hamilton  &  Healy  1969;  Ryall  &  Savage  1969)  as  well  as 
studies  of  faulting  and  fracturing  caused  by  the  explosions 
(Buckman  1969;  McKeown  &  Dickey  1969)  have  been  published. 

Figure  10  shows  the  relationships  of  these  two  sets  of  data  with 
the  source  mechanisms  determined  for  the  two  events.  In  the 
case  of  Benham,  the  correlation  between  the  orientation  of  the 
double-couple  and  the  trend  of  the  aftershock  pattern  as  well 
as  the  orientation  of  faults  and  fractures  is  quite  good.  Both 
vertical  and  right-lateral  movements  on  the  nearby  Boxcar  fault 
were  initiated  by  Benham. 

The  aftershocks  of  Benham  were  concentrated  in  two  patterns 
(Hamilton  &  Healy  1969).  Those  to  the  northwest  of  ground  zero 
follow  a  northeasterly  trend  while  those  to  the  west  and  south¬ 
west  trend  north-south.  Fault  plane  solutions  from  P  waves 
of  events  in  the  first  group  indicate  dip-slip  fault  movement. 
Solutions  for  the  second  group  indicate  right-lateral  strike- 
slip  motion  in  a  northerly  direction.  These  two  patterns  were 
further  confirmed  in  a  study  of  smaller  Benham  aftershocks 
(Stauder  1971) .  The  solutions  for  the  second  group  of  after¬ 
shocks  are  similar  to  the  orientation  of  the  double-couple  com¬ 
ponent  for  Benham  (and  other  events)  as  determined  by  this  study. 
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The  Benham  and  Boxcar  aftershocks  when  considered  together 
(see  Fig.  10)  indicate  the  continuity  of  the  distribution  pat¬ 
tern.  Northeast-trending  Benham  aftershocks  overlap  with  those 
of  Boxcar . 

II.  EFFECT  OF  TECTONIC  STRAIN  RELEASE  ON  MAGNITUDES  AND 
DISCRIMINANTS 

With  explosions  such  as  Pile  Driver,  Hardhat  and  Greeley 

where  the  surface  wave  energy  due  to  strain  release  is  greater 

than  that  of  the  explosion,  one  might  expect  the  M  values 

b 

to  increase.  The  expected  theoretical  increase  AM  deter- 

b 

mined  from  the  Gutenberg- Richter  energy-magnitude  relationship 
for  different  F  values  is  listed  in  Table  2.  These  were  com¬ 
puted  using 

AMS  =  ^  log  [l  +  Etect/Eexp]  (5) 

For  explosions  such  as  Pile  Driver  and  Greeley  (with  F 
values  of  3.2  and  1.6,  respectively)  one  would  expect  surface 
wave  magnitudes  to  be  greater  than  those  of  corresponding  explo¬ 
sions  with  small  F  values.  One  way  of  investigating  this 
point  is  to  look  at  Mg  vs  m^  values  for  a  series  of  explosions 
with  different  F  values.  In  Fig.  11,  Mg  -  explosion  data 


-85- 


from  three  different  authors  are  plotted.  Neither  Pile  Driver 
nor  Greeley  is  distinguishable  from  the  general  explosion  popu¬ 
lations.  The  reason  for  this  seems  to  be  that  Rayleigh  waves 
due  to  the  explosion  component  and  double-couple  component  of 
the  source  reinforce  and  cancel  out  in  alternate  azimuthal  qua¬ 
drants.  Thus  when  azimuthally  averaged,  the  net  effect  on  mag¬ 
nitudes  is  small. 

To  illustrate  this  point  let  us  compare  the  Pile  Driver 
radiation  pattern  and  azimuthal  range  of  the  Canadian  network 
stations  used  by  Basham  (1969)  in  his  magnitude  determinations 
in  Fig.  12.  With  our  source  model,  Canadian  stations  can  be 
expected  to  record,  on  the  average,  normal  Rayleigh  wave  ampli¬ 
tudes  and  thus  yield  a  normal  vs  m^  value,  as  is  indeed 
observed. 

Thus  to  insure  a  representative  magnitude,  it  is  important 
to  use  data  from  several  azimuths.  Furthermore,  while  establish¬ 
ing  station  corrections  for  magnitude  determinations,  it  is 
important  to  use  data  from  explosions  with  low  F  values  or 
from  those  at  different  sites. 


III.  SOURCE-TIME  FUNCTION 

All  observations  that  have  been  made  on  the  spectra  of 
surface  waves  from  explosions  and  earthquakes  of  comparable 


fiwssi-aiCT&tfU 


-86- 


magnitudes  indicate  that  earthquake  spectra  are  richer  in  long- 
period  components  compared  to  explosions.  This  could  be  due 
to  a)  source  depth,  b)  source  volume  and  c)  source- time  function. 

Tsai  &  Aki  (1971)  favor  the  source-depth  effect  to  ex¬ 
plain  the  differences  in  spectra  while  assuming  the  source 
function  for  explosions  and  earthquakes  to  be  the  same  (step 
function) .  The  spectra  of  deeper  earthquakes  would  be  shifted 
toward  lower  frequencies.  This  effect  cannot  be  neglected  but 
it  does  not  explain  two  sets  of  data:  1)  Some  aftershocks  of 
the  Benham  and  Jorum  events  were  accurately  located.  While 
focal  depths  were  lesr  than  5  km  (some  as  shallow  as  3.5  km) 
their  spectral  characteristics  were  similar  to  those  of  earth¬ 
quakes  rather  than  explosions  (Molnar  et  al  1969;  Savino  et  al 
1971) .  Also,  spectra  of  collapse  events  are  shifted  toward 
lower  frequencies  relative  to  explosions.  2)  Spectra  of  Love 
waves  generated  by  the  explosions  are  similar  to  those  generated 
by  earthquakes  (Savino  et  al  1971;  ToksdJz  et  al  1971).  In  fact 
the  Love  wave  source-time  function  determined  from  amplitude 
equalization  is  close  to  a  step  function. 

For  an  explosion,  a  decaying  source-pressure  function, 
with  or  without  a  residual  d-c  pressure,  could  explain  the 
observed  spectra.  In  our  earlier  studies  (Toksdz  et  al  1969) 
we  adopted  a  pressure  function  given  by 

p(t)  =  PQte-nt 


(6) 
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where  pQ  is  a  constant  and  ri  is  a  time  constant  dependent 
on  yield  in  a  given  medium.  For  Bilby,  n  =1.5  sec  ^  fits 
the  observed  spectra  at  three  stations  as  ir.  shown  in  Fig.  13. 
Discrepancies  at  periods  longer  than  T  =  30  sec  are  attributed 
to  the  effects  of  lateral  heterogeneities  and  low  signal-to- 
noise  ratio. 

In  a  given  medium,  the  source  pulse  becomes  broader  with 
increasing  yield.  This  could  be  explained  by  the  larger  non¬ 
linear  zone  and  hence  greater  attenuation  of  higher  frequencies 
in  this  zone  (Andrews  &  Shlien  1972).  Observed  ground  displace¬ 
ment  spectra  of  Rayleigh  waves  from  several  NTS  explosions  as 
recorded  at  LASA  is  shown  in  Fig.  14  along  with  three  theoreti¬ 
cal  spectra  for  two  source  types.  Curve  1  is  based  on  a  step 
function  pressure  variation.  Curves  2  and  3  are  based  on  the 
decaying  pressure  pulse  with  rj  =  1  and  r}  =0.6  sec  \  respec¬ 
tively.  It  is  clear  that  for  the  five  lower  explosions  with  in¬ 
termediate  yields,  the  decaying  source  function  with  r)  =1.0 
fits  the  spectra,  while  for  the  three  larger  events,  n  <  1  is 
necessary.  The  specra  of  the  step  function  and  the  decaying 
pulse  [Eq.  (6)]  type  source  functions  are  given  by: 


|S(W)  |  =  ~  T/2-n 


step 


(7) 


|  S  (w)  | 


decaying  pulse 


-88- 


where  w  is  the  angular  frequency  and  T  is  the  period.  At 
very  low  frequencies,  the  spectrum  of  the  decaying  pulse 
approaches  a  constant,  while  that  of  the  step  increases  with 
increasing  period.  At  higher  frequencies  the  pulse  spectrum 
falls  off  proportionally  to  w  ,  while  the  step  falls  off  as 
1/tu. 

The  success  of  the  spectral  ratio  discriminant  (ratio  of 
amplitude  at  T  =  20  sec  to  T  —  40-60  sec)  can  also  be  ex¬ 
plained  by  the  difference  in  source-time  functions  of  the  ex¬ 
plosions  and  earthquakes  as  given  in  Eq.  (7). 
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CONCLUSIONS 

Nuclear  explosions  detonated  in  media  with  tectonic  stresses 

interact  with  the  ambient  stress  field  and  release  some  of  the 

strain  energy.  This  phenomenon  complicates  the  Rayleigh  wave 

radiation  patterns  and  generates  Love  waves.  It  could  also 

affect  the  -  m.  discriminants.  Some  conclusions  that  can 
d  b 

be  derived  from  this  study  are: 

1.  The  strain  energy  release  from  a  particular  explosion 
is  dependent  on  the  rock  type  and  the  state  of  stress  in  the 
medium.  Generally,  explosions  in  harder  rocks  rv lease  greater 
amounts  of  strain  energy.  This  probably  is  related  to  the 
energy-storing  capability  of  these  rocks.  No  significant  Love 
waves  have  been  observed  from  explosions  in  loose  alluvium  or 
salt  domes. 

2.  In  a  given  rock  type  tectonic  strain  energy  release 
varies.  As  was  shown  by  laboratory  experiments  (Toksflz  el  al 
1971),  this  is  related  to  the  ambient  stress  level. 

3.  The  strength  of  the  tectonic  component  of  surface 
waves  relative  to  the  explosion  itself  increases  with  period. 

4.  Tectonic  strain  energy  released  as  surface  waves  is 
generally  less  than  the  surface  wave  energy  due  to  the  explosion 
alone.  In  three  cases  (Pile  Driver,  Hardhat,  and  Greeley)  the 
tectonic  component  of  the  surface  waves  exceeded  those  due  to 
the  explosion  and,  therefore,  dominated  the  radiation  pattern. 


5.  The  tectonic  strain  release  characteristics  of  explo¬ 
sions  did  not  have  significant  effect  on  the  Mg  -  m^  discri¬ 
minant.  This  may  be  due  to  two  factors:  a)  The  mean  Rayleigh 
wave  amplitude  when  averaged  over  all  azimuths,  does  not  change 


as  a  result  of  the  tectonic  strain  release  component.  Explosion 
and  tectonic  components  of  the  source  interfere  constructively 
and  destructively  in  alternating  quadrants,  b)  Generally,  for 
an  explosion  of  given  yield,  the  body  wave  amplitudes  are  larger 
when  the  source  is  in  harder  media.  Azimuthal  averaging  is  most 
important  for  Mg  in  Mg  -  m^  discrimination. 

6.  A  decaying  pulse-type  source-time  function  can  explain 
surface  wave  amplitude  spectra  and  long-period/ short-period 
spectral  ratio.  The  shape  of  the  pulse  depends  on  medium  pro¬ 
perties  and  size  of  the  explosion. 
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TABLE  1 . 

TECTONIC  STRAIN  RELEASE  CHARACTERISTICS  OF  UNDERGROUND  NUCLEAR 

EXPLOSIONS 


U.  S 


EVENT 

REGION 

Pile  Driver 

6-2-66 

Yucca  Flat 
(n.  end) 

Hardhat 

2-15-52 

Yucca  Flat 
(n.  end) 

Shoal 

10-26-63 

Fallon, 

Nevada 

Greeley 

12-20-66 

Pahute  Mesa 

Benham 

12-19-68 

Pahute  Mesa 

Chartreuse 

5-6-66 

Pahute  Mesa 

Duryea 

4-14-66 

Pahute  Mesa 

Half  Beak 

6-30-66 

Pahute  Mesa 

Boxcar 

4-26-68 

Pahute  Mesa 

Corduroy 

12-3-65 

Yucca  Flat 

Rulison 

9-10-69 

Grand  Valley, 
Colorado 

Faultless 

1-19-68 

Central 

Nevada 

Cup 

3-26-65 

Yucca  Flat 

Bilby 

9-13-63 

Yucca  Flat 

Tan 

6-3-66 

Yucca  Flat 

Bronze 

7-23-65 

Yucca  Flat 

Buff 

12-16-65 

Yucca  Flat 

Haymaker 

6-27-62 

Yucca  Flat 

Sedan 

7-6-62 

Yucca  Flat 

Salmon 

10-22-64 

Hattiesburg, 

Mississippi 

Gnome 

12-10-61 

Carlsbad, 

New  Mexico 

Milrow 

10-2-69 

Amchitka 

Cannikin 

11-5-71 

Amchitka 

Explosions 


MEDIUM 

D.C. 

STRENGTH 

(F) 

FAULT 

AZI. 

ENERGY  RATIO 

E,  ,/E 

tect  exp 

Granite 

3.20 

340° 

13.65 

Granite 

3.00 

330° 

12.00 

Granite 

.90 

346° 

1.05 

Zeolite 

1.60 

355° 

3.41 

Tuff 

Zeolite 

.85 

345° 

.96 

Tuff 

Rhyolite 

,90 

353° 

1.05 

Rhyolite 

.75 

355° 

.75 

Rhyolite 

.67 

345° 

.60 

Rhyolite 

.59 

346° 

.46 

Quartzite 

.72 

347° 

.69 

Ss.  and 

.60 

335° 

.48 

Shale 

Sat.  Tuff 

.50 

344° 

.33 

Tuff 

.55 

200° 

.40 

Tuff 

.47 

340° 

.29 

Tuff 

.39 

347° 

.20 

Tuff 

.33 

185° 

.15 

Tuff 

.31 

208° 

.13 

Alluvium 

.33 

340° 

.14 

Alluvium 

0 

- 

0 

Salt 

0 

- 

0 

Salt 

0 

- 

0 

Andesite 

<  .6 

- 

<.48 

Andesite 

.36 

- 

.17 

\ 
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TABLE  1 .  Continued 
Presumed  Soviet  Explosions 


DATE 


EVENT 

(m-d-y) 

REGION 

U.S.S.R. 

10-27-66 

Novaya  Zemlya 

U.S.S.R. 

10-21-67 

Novaya  Zemlya 

U.S.S.R. 

2-26-67 

E.  Kazakh 

U.S.S.R. 

3-20-67 

E.  Kazakh 

U.S.S.R. 

2-13-66 

E .  Kazakh 

MEDIUM 

D.C. 

STRENGTH 

(F) 

FAULT 

AZI. 

ENERGY  RATIO 
Etect^Eexp 

-- 

0.90 

5° 

1.05 

-- 

0.71 

43° 

0.67 

-- 

0.85 

6° 

0.96 

0.81 

155° 

0.87 

_  _ 

0.67 

101° 

0.60 

EFFECT  OF 


F 

Double-Couple 

Part 

0.3 

0.5 

0.7 

1.0 

1.5 

2.0 

3.0 

4.0 


TABLE  2 


STRAIN  RELEASE  ON  MAGNITUDES 


E.  . /E 
tect'  exp 


AMS 

Maximum  increase  in  magnitude 
due  to  tectonic  component 


0.12 

0.33 

0.65 

1.33 
3.00 

5.33 
12.00 
21.33 


0.07 

0.08 

0.14 

0.24 

0.40 

0.53 

1.08 

1.23 


Figure  1. 


Figure  2. 


Figure  3. 


Figure  4. 


Figure  5. 
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Figure  Captions 

Long-period  filtered  seismograms  from  the 
Pile  Driver,  Tan,  and  Greeley  explosions 
recorded  at  Weston,  Massachusetts. 

Stations  recording  surface  waves  from  explosions 
in  the  continental  United  states. 

Filtered  (upper  three  traces)  and  unfiltered 
(lower  three  traces)  long-period  components  from 
the  Boxcar  explosion  at  Mould  Bay,  Canada  (MBC) 
pbtted  from  digital  data  at  the  same  gain.  The 
two  filtered  horizontal  components  have  been 
rotated  to  radial  and  tangential  directions  to 
separate  Love  and  Rayleigh  waves . 

Contour  plots  of  the  deviations  of  combinations 
of  part  double-couple  and  fault  plane  azimuth 
from  experimental  Love  to  Rayleigh  wave  ratios 
for  three  period  ranges:  a)  9  to  15  sec., 
b)  15  to  22  sec. ,  and  c)  22  to  30  sec. ,  for  the 
Faultless  event. 

The  period  dependence  of  F,  the  part  double¬ 
couple,  for  five  nuclear  explosions. 


Figure  6. 


Figure  7. 
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Observed  and  theoretical  Love/Rayleigh  wave 
anplitude  radiation  patterns  for  three  explosions 
with  different  proportions  of  strain  release 
(F)  . 

Raleigh  wave  radiation  patterns  for  explosive 
component  (circles) ,  strain  release  component 
(broken  curves) ,  and  composite  sources 
representing  the  theoretical  models  (solid 
curves,  where  for  Pile  Driver:  F  =  3.2  and  ©  =  340°; 
and  for  Greeley:  F  =  1.6  and  0  =  355°) .  The 
polarity  (phase)  of  Rayleigh  waves  in  different 
quadrants  of  the  radiation  pattern  is  indicated 
by  (+)  and  (-) .  For  Pile  Driver  some  sample 
seismograms  from  four  stations  are  compared  to 
those  of  Tan.  Amplitude  factors  are  arbitrary. 

Pile  Driver  traces  are  below  those  of  Tan. 

Dashed  line  traces  have  their  polarities  reversed. 
Note  the  perfect  match  of  wave  shapes,  with 
polarities  reversed  at  JP-AT  and  BOZ  as  pre¬ 
dicted  by  the  model.  For  Greeley  sample  seismograms 
are  compared  to  those  of  Half  Beak  (COL  and  OXF) 
and  Boxcar  (LASA) .  The  lower  traces  are  from 
Greeley.  Polarity  reversal  is  very  clear  at 
LASA.  MN-NV,  Mina,  Nev. ; JP-AT,  Jasper,  Alta.; 

BOZ,  Bozman,  Mont.;  AX2AL ;  Alexander  City,  Ala.; 


Figure 


Figure 


Figure 
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TUC,  Tucson,  Ariz.;  COL,  College,  Alaska;  LASA, 
Large  Aperature  Seismic  Array,  Mont.;  OXF ,  Oxford, 
Miss.  Data  symbols:  A LRSM  stations;  A,  LRSM 
stations  where  reversed  polarity  was  observed;  ■ 
WWSSN  stations;  •,  CSSN  stations. 

8.  Map  of  Yucca  Flat  showing  the  fault  plane  orienta¬ 
tion  of  the  explosions  studied,  explosion-produced 
fractures,  and  major  natural  faults  in  the  area 
(geology  after  Fernald  et  al  1968) . 

9.  Map  of  Pahute  Mesa  showing  the  fault  plane 
orientations  of  the  explosions  studied,  explosion- 
produced  fractures,  and  natural  faults  in  the 
area  (faults  after  Orkild  et  al  1969) . 

10.  Source  mechanisms  determined  for  the  Boxcar  and 
Benham  events  compared  to  natural  faults  and 
explosion-produced  fracturing,  and  aftershock 
distribution.  Aftershock  distribution  for  Benham 
(Hamilton  &  Healy  1969) , is  that  of  the  larger 
aftershocks  (magnitude  3  to  4.2).  Aftershocks 
from  Boxcar  (Ryall  &  Savage  1969)  are  those 
occurring  in  the  first  two  days  after  the  event. 
The  circles  are  of  radius  5  km  centered  on  the 
ground  zero  of  the  explosion. 
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Figure  11.  Mg  -  i^b  for  13  American  underground  nuclear 

explosions  from  the  data  of  Liebermann  &  Pomeroy 
(1969)  ,  Basham  (1969),  and  Wagner  (1970). 

Figure  12.  Rayleigh  wave  radiation  pattern  for  Pile  Driver 
with  the  azimuthal  range  of  Canadian  stations 
from  the  event. 

Figure  13.  Amplitude  spectra  of  the  source-time  function 

for  Bilby  (n  *1.5  sec-1)  at  three  stations  after 
correction  for  propagation  effects.  Circles  in¬ 
dicate  data  uncorrected  for  attenuation;  triangles 
indicate  the  data  corrected  taking  Q  =  100. 

Figure  14.  Vertical  component  amplitude  spectra,  corrected 
for  instrument  response,  from  eight  explosions 
recorded  at  LASA.  Curves  1,  2,  and  3  are  computed 
Rayleigh  wave  spectra  (observed  explosion  spectra 
and  Curves  1  and  2  from  Lande  &  Filson  1970) . 
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4 .  EARTH  STRUCTURE 

4.1  Seismic-Wave  Attenuation  and  Partial  Melting  in 

the  Upper  Mantle  of  North  America  by  Sean  C .  Solomon 

(Abstract) 

A  model  of  Q*”1  based  on  Walsh's  theory  for  attenuation 
in  partially  melted  rock  is  proposed  for  the  upper  mantle 
of  western  North  America.  The  asthenosphere  (or  low-Q  zone), 
in  which  attenuation  is  attributed  to  a  superposition  of 
thermally  activated  relaxation  processes,  is  300  km  thick  in 
the  irodel  and  must  be  vertically  inhomogeneous.  The  litho¬ 
sphere  (or  high-Q  lid)  is  60  km  thick.  The  model  is  con¬ 
sistent  with  a  wide  assortment  of  attenuation  and  velocity 
da*a,  in  all  spanning  approximately  three  decades  in  frequency. 
Both  the  shear  modulus  and  Q-1  for  shear  waves  are  frequency 
dependent  in  the  asthenosphere.  The  viscosity  and  volume 
concentration  of  melt  can  be  estimated  from  the  relaxation 
parameters.  Further,  if  lateral  changes  in  these  parameters, 
determined  from  the  differential  attenuation  and  travel-time 
delays  of  P  and  S  waves,  are  attributed  to  horizontal 
temperature  variations,  then  these  temperature  differences  in 
the  asthenosphere  need  be  no  larger  than  100  to  200°K  at  a 
depth  of  several  hundred  kilometers.  The  highest  temperatures 
determined  in  this  fashion  are  beneath  the  Rocky  Mountain 
and  Pacific  border  regions. 
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4.2  On  Q  and  Seismic  Discrimination  by  Sean  C.  Solomon 
Summary 

The  vagaries  of  seismic  attenuation,  or  Q,  in  the  earth 
render  uncertain  the  estimates  of  many  seismic  source  para¬ 
meters  made  by  the  unwary.  In  particular,  the  discrimination 
of  waves  produced  by  earthquakes  from  those  of  explosicns,  if 
the  waves  have  traveled  dissimilar  paths,  or  the  characteri 
zation  of  a  seismic  source  from  the  seismic-wave  amplitude 
spectrum  require  concern  for  the  spatial  and  frequency  depen¬ 
dence  of  Q.  Surface-wave  amplitudes  are  most  affected  by 
attenuation  at  periods  long  enough  (30  sec  and  greater)  so 
that  much  of  the  energy  propagates  in  the  low-Q  as theno sphere, 
and  at  periods  short  enough  (less  than  15  sec)  so  that  near¬ 
surface  structure  is  important.  Thus  the  standard  surface- 
wave  magnitude  is  largely  insensitive  to  anelastic  losses, 
though  for  the  surface-wave  spectral  shape  to  be  diagnostic 
of  source  depth,  such  losses  must  be  considered.  The  ampli¬ 
tudes  of  body  waves  that  have  penetrated  the  low-Q  astheno- 
sphere  reflect  the  extreme  variations  of  Q  within  that  zone. 
Seismic  discrimination  criteria  based  on  the  body-wave 
magnitude  should  be  restricted  to  comparison  of  events  from 
the  same  tectonic  province.  Further,  estimates  of  such  source 
properties  as  source  dimension  and  fractional  stress  drop 
depend  critically  on  properly  accounting  for  the  frequency 
dependence  of  Q. 
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1 .  Introduction 


The  problem  of  discriminating  shallow  earthquakes  from 
underground  nuclear  explosions  may  be  visualized  as  an  attempt 
to  distinguish  between  the  signals  from  two  dissimilar,  broad¬ 
band  sources  after  those  signals  have  passed  through  a  compli¬ 
cated  narrow-band  filter,  the  earth.  That  the  characteristics 
of  this  terrestrial  filter  depend  upon  the  precise  location 
of  source  and  receiver  adds  further  difficulties.  A  major 
component  of  the  filter  is  seismic-wave  attenuation,  most  con¬ 
veniently  represented  by  the  quantity  Q,  the  dimensionless 
quality  factor.  The  distribution  of  Q  in  the  earth  determines 
both  the  absolute  amplitudes  of  the  various  types  of  seismic 
waves  and  the  shape  of  their  amplitude  spectra.  In  this  paper, 
the  measured  properties  of  seismic  attenuation  in  the  earth  are 
shown  to  have  potentially  important  effects  on  estimates  of 
seismic  source  parameters  and,  in  particular,  on  some  of  the 
techniques  currently  used  to  discriminate  between  earthquakes 
and  explosions. 

It  should  perhaps  be  stated  at  the  outset  that  the 
remarks  below  have  no  special  bearing  on  the  seismic  discrim¬ 
ination  problem  if  explosion  and  earthquake  source  are  in 
close  proximity  and  if  signals  from  both  are  observed  by  the 


\ 


same  receiver.  Only  if  signals  that  have  traveled  different 
paths  or  that  have  been  recorded  over  different  frequency 
bands  are  to  be  compared  or  contrasted  does  careful  note 
need  be  made  of  corrections  for  attenuation  and  other 
phenomena  of  propagation. 


2 .  Some  Properties  of  Q 

There  are  three  features  of  the  distribution  of  Q  in 
th*3  earth  that  are  relevant  to  the  discussion  of  seismic 
discrimination : 

(1)  Depth  dependence.  Q,  considered  as  a  material 
property,  is  a  strong  function  of  depth  in  the  earth.  Studies 
bohh  of  the  attenuation  of  surface  waves  and  free  oscillations 
as  a  function  of  period  (Anderson  et  al.  1965;  Tsai  &  Aki  1969) 
and  of  body-wave  amplitudes  as  a  function  of  distance  (Wadati 
&  Hirono  1956;  Archambeau  et  al.  1969)  indicate  a  zone  in  the 
upper  mantle  within  which  Q  is  significantly  lower  than  in  the 
crust  or  the  remainder  of  the  mantle.  This  low-Q  zone  is  iden¬ 
tified  with  the  asthenosphere  in  the  new  global  tectonics 
(Isacks  et  al.  1968) .  Typically  Q  in  the  asthenosphere  is  an 
order  of  magnitude  lower  than  in  the  overlying  lithosphere. 
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(2)  Lateral  variation.  Q  is  also  region-dependent.  The 
most  dramatic  lateral  variations  occur  at  island  arcs,  where 
the  high-Q  lithosphere  is  downthrusting  into  the  mantle  (Oliver 
&  Isacks  1967)  ,  and  at  other  regions  where  the  l  ithosphere  is 
abnormally  thin  or  discontinuous  (Molnar  &  Oliver  1969) .  Q 
also  shows  lateral  changes  within  the  lithosphere  and  astheno- 
sphere.  The  most  pronounced  effects  on  seismic-wave  amplitudes 
recorded  at  different  stations  are  produced  by  the  Q  variations 
within  the  low-Q  asthenosphere  (Solomon  &  Toksftz  1970) .  These 
variations  are  most  simply  explained  by  relatively  modest 
horizontal  temperature  gradients  (Solomon  1972) . 

(3)  Frequency  dependence.  At  least  within  the  astheno- 
sphere,  Q  is  a  function  of  frequency.  That  Q  may  increase 
with  frequency  for  short-period  P  waves  in  the  mantle  is 
suggested  by  the  high  values  of  Qp  (greater  than  1000)  deter¬ 
mined  for  the  frequency  range  1-10  Hz  (Asada  &  Takano  1963; 
Sumner  1967;  Filson  1970)  contrasted  with  values  an  order  of 
magnitude  lower  for  frequencies  in  the  range  .01  to  1  Hz 
(Teng  1968;  Hirasawa  &  Takano  1966;  Kanamori  1967) .  Further, 

Q  approximately  proportional  to  frequency  has  been  observed 
for  P  and  S  waves  over  various  paths  and  over  assorted  fre¬ 
quency  bands  within  the  frequency  range  .3  to  20  Hz  (Frantti 
1965;  Kurita  1968;  Fedotov  &  Boldyrev  1969;  Archambeau  et  al. 
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1969;  Mendoza  1970).  The  suggestion  has  also  been  made  that 
the  attenuation  of  surface  waves  and  free  oscillations  is 
best  fit  by  a  frequency-dependent  Q  in  the  mantle  (Tsai  & 

Aki  1969;  Jackson  1972). 

The  best  test  of  the  frequency  dependence  of  Q  comes 
from  a  consideration  of  data  taken  over  a  broad  frequency 
band.  A  set  of  measurements  of  elastic  and  anelastic 
properties  of  the  upper  mantle  of  western  U.S.,  in  all 
spanning  a  factor  of  1000  in  frequency,  could  not  all  be 
satisfied  by  a  Q  model  independent  of  frequency  (Solomon  1972) . 
On  the  other  hand,  a  relatively  simple,  frequency-dependent 
Q  model,  based  on  the  expected  behaviour  of  attenuation  in 
partially  melted  rock  (Walsh  1969) ,  could  be  made  to  satisfy 
the  entire  set  of  observations. 


3.  Surface  waves 


The  above  properties  of  seismic  attenuation  in  the  earth 
affect  surface  waves  and  body  waves  differently,  so  these 
two  modes  o t  energy  propagation  are  considered  separately 
below.  For  our  purposes,  the  important  distinctions  between 
surface  and  body  waves  are  the  typically  fewer  number  of 


/ 
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wavelengths  per  total  propagation  path  and  the  generally 
narrower  width  of  the  relevant  frequency  band  of  the  former 
relative  to  the  latter.  These  distinctions  imply  that  the 
amplitudes  of  Love  and  Rayleigh  waves  from  a  given  seismic 
event  are  less  controlled  by  anelastic  losses  than  are  those 
of  P  and  S  waves,  and  that  surface-wave  amplitudes  are 
generally  less  sensitive  than  those  of  body  waves  to  the  fre¬ 
quency  dependence  of  Q. 

Principally  because  Q  in  the  earth  is  a  strong  function 
of  depth,  the  attenuation  of  surface  waves  is  period -dependent . 
The  contribution  of  any  depth  interval  in  the  earth  to  the 
attenuation  of  a  surface  wave  is  proportional  to  the  frac¬ 
tional  potential  energy  stored  per  cycle  within  that  interval 
(Anderson  &  Archambeau  1964) .  In  Figure  1  is  shown  the  frac¬ 
tional  potential  energy  per  unit  depth  stored  per  cycle  for 
Rayleigh  waves  of  10,  20  and  40  sec  periods.  The  boundary 
between  high-Q  lithosphere  and  low-Q  asthenosphere  is  indicated 
schematically  at  a  depth  of  60  km.  Discontinuities  in  the 
energy  as  a  function  of  depth  are  due  to  abrupt  changes  in 
the  elastic  properties  of  the  particular  velocity-density 
model  used,  the  shield  model  of  Harkrider  &  Anderson  (1966) . 

The  important  point  of  Figure  1  is  that  Rayleigh  waves  of  10- 
and  20-sec  period  store  a  relatively  small  fraction  (.002  and 
5  percent,  respectively)  of  their  strain  energy  in  the  low-Q 
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asthenosphere,  while  40-sec  Rayleigh  waves  store  over  half 
their  potential  energy  per  cycle  below  60  Ion  depth.  Thus 
in  regions  of  the  earth  where  Q  in  the  asthenosphere  is 
low,  the  attenuation  of  surface  waves  of  40  sec  period  should 
be  more  pronounced  than  at  shorter  periods .  At  periods  less 
than  about  10  to  15  sec,  surface-wave  attenuation  is  very 
sensitive  to  near-surface  values  of  Q  (see  Figure  1)  which 
can  also  be  path-dependent.  For  oceanic  or  continental  paths 
crossing  sediment  basins,  for  instance,  surface-wave  attenu¬ 
ation  at  these  periods  might  be  expected  to  be  severe. 

In  Figures  2  and  3  are  given  the  measured  Love-wave 
and  Rayleigh-wave  attenuation  (QL“1  and  QR_1)  for  east-central 
and  western  United  States  and  for  predominantly  ocean  paths. 
The  values  for  western  U.S.  were  reported  earlier  (Solomon 
1972) ;  those  for  east-central  U.S.  were  determined  using  stan¬ 
dard  two-station  techniques  from  surface  waves  which  had 
traveled  between  WWSSN  stations  at  Rapid  City,  South  Dakota 
(RCD)  and  Atlanta,  Georgia  (ATL) .  Details  of  both  sets  of 
measurements  may  be  found  in  Solomon  (1971) . 

At  periods  less  than  25  to  30  seconds,  Q  and  Q  exceed 

R 

200  and  are  not  markedly  path-dependent.  This  reflects  the 
uniformly  high  Q  of  the  lithosphere.  At  longer  periods,  both 
Qr  1  and  Ql  1  generally  increase  with  increasing  period,  due 
to  the  relatively  larger  fraction  of  energy  stored  in  the 
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asthenosphere  at  these  periods.  Further,  since  lateral  vari¬ 
ations  in  the  anelastic  properties  of  the  earth  are  pronounced 
within  the  asthenosphere,  the  regional  differences  in  surface- 
wave  attenuation  are  greatest  in  the  period  range  30  to  about 
100  sec.  Q.”1  and  q”1  for  the  tectonically  active  western 
U.S.  are  consistently  larger  in  this  period  range  than  for 
the  more  stable  regions  to  the  east. 

These  findings  have  considerable  significance  for  the 
interpretation  of  surface-wave  magnitudes  and  spectra.  In 
Figure  4  is  shown  the  vertical  displacement  amplitude  spec¬ 
trum  of  a  Rayleigh  wave  from  a  nuclear  explosion  (Tsai  &  Aki 
1971) .  Also  shown  are  the  Rayleigh-wave  spectra  predicted 
for  3000-km  propagation  paths  through  "western  U.S."  and 
"eastern  U.S."  (A  flat  earth  has  been  assumed,  so  that  geo¬ 
metrical  spreading  is  ignored.)  At  periods  less  than  25  sec, 
the  attenuated  spectra  are  quite  similar.  In  particular,  the 
surface-wave  magnitude,  if  computed  from  the  amplitude  of 
20-sec  waves  (Gutenberg  1S45) ,  would  be  identical  for  western 
and  eastern  U.S.  paths. 

At  periods  longer  than  25  sec,  the  two  attenuated  spectra 
diverge.  Thus  some  concern  for  propagation  path  must  be 
exercised  if  one  bases  a  discrimination  criterion  on  the 
amplitudes  of  40  sec  surface  waves  (Molnar  et  al.1970;  Savino 
et  al.  1971).  Still,  differences  in  surface-wave  "magnitude" 
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would  amount  to  less  than  .1  at  periods  of  40  to  50  sec, 
and  less  than  .3  at  periods  near  30  sec  for  the  two  3000-km 
paths. 

Anothor  cautionary  note  should  be  raised.  Anelastic 
losses  do  affect  the  shape  of  surface-wave  amplitude  spectra 
(Figure  4) .  Thus  discrimination  criteria  based  on  the  sensi¬ 
tivity  of  surface— wave  spectral  shape  to  depth  of  focus  (Tsai 
&  Aki  1971)  require  accurate  knowledge  of  Q.  The  attenuated 
spectra  of  Figure  4,  if  "corrected"  for  attenuation  by  assuming 
Q  is  proportional  to  frequency  (Tsai  &  Aki  1970,  1971),  would 
appear  to  be  from  an  event  as  deep  as  5  km  (see  Figures  5,  6 
and  7  of  Tsai  &  Aki  1971)  .  Such  an  event  might  be  incorrectly 
classified  as  an  earthquake  if  good  azimuthal  coverage  were 
not  possible  or  if  other  discrimination  criteria  were  not 
available. 


4 .  Body  waves 

The  amplitudes  of  P  and  S  waves  depend  criticall^  on 
both  the  spatial  and  frequency  dependence  of  Q  in  the  earth. 
Because  Q  in  the  lithosphere  is  high,  the  attenuation  correction 
for  P  and  S  waves  recorded  within  a  few  hundred  km  of  the  source 
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will  be  small  and  insensitive  to  regional  variations.  For 
source-receiver  distances  greater  than  about  1000  km, 
direct  P  and  S  waves  must  propagate  through  the  asthenosphere 
or  low-velocity  zone,  in  which  the  waves  are  strongly  damped 
and  within  which  lateral  variations  in  Q  can  noticeably  alter 
the  relative  amplitudes  between  different  stations. 

One  measure  of  the  lateral  variation  of  Q  in  the  astheno¬ 
sphere  of  North  America  is  the  differential  attenuation  of 
long-period  body  waves  (Solomon  &  Toksflz  1970)  ,  mapped  for 
P  waves  in  Figure  5.  The  differential  attenuation  fit*  at  a 
seismometer  station  is  related  to  the  difference  between  the 
value  of  Q-^  in  the  upper  mantle  beneath  the  station  and 
the  average  Q-^  for  a  spherically  symmetric  earth: 

fit*  =  tfJiq”1  -  <fr] v_1ds 
S 

where  ds  is  an  increment  of  distance  along  the  ray  path  S 
and  v  is  the  wave  velocity.  In  practice,  the  difference 
between  fit*  at  two  stations  is  given  by  the  negative  of 
the  slope  of  the  body-wave  spectral  ratio,  corrected  for 
the  different  distances  from  the  source  to  the  two  stations, 
versus  frequency.  Sources  of  body  waves  used  to  determine 
values  of  fi  t*  in  Figure  5  were  deep  earthquakes  at  distances 
greater  than  40  degrees,  so  that  variations  in  Q  along  near-source 
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segments  of  paths  were  minimized. 

In  Figure  5  may  be  seen  a  large  area  of  above-average 
attenuation  (positive  St*)  of  P  waves  between  the  Rocky 
Mountain  front  and  the  Sierra  Nevada-Cascade  ranges.  One 
consequence  of  Figure  5,  if  variations  in  upper  mantle  Q  ^ 
implied  by  differences  in  St*  may  be  at  least  qualitatively 
extrapolated  to  higher  frequencies,  is  that  the  magnitudes 
of  P  waves  recorded  from  teleseismic  events  at  stations  in 
the  Basin  and  Range  and  Rocky  Mountain  provinces  should  be 
smaller  than  those  determined  at  eastern  U.S.  stations.  This 
is  generally  observed  (Cleary  1967;  Evernden  &  Clark  1970). 

A  second  consequence  is  that  the  body-wave  magnitude  m^  of 
earthquakes  in  the  Basin  and  Range  province  should  be  less 
than  for  earthquakes  of  comparable  surface  wave  magnitude 
M  in  the  neighboring  California  Border  and  Great  Plains 

b 

provinces,  providing  that  body  waves  are  recorded  at  tele- 
seismic  distances,  that  surface  waves  of  period  20  sec  or 
less  are  used  to  determine  Mg,  and  that  depth  of  focus 
and  other  source  parameters  are  similar  for  earthquakes  in 
the  two  regions . 

As  a  test  of  this  latter  hypothesis,  we  reproduce 
in  Figure  6  Basham's  (1969)  measurements  at  Canadian  stations 
of  and  Mg  for  earthquakes  and  explosions  in  southwestern 

North  America.  Unlike  Basham's  original  figure,  the  earthquake 
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population  is  segregated  into  two  groups:  (i)  those  that 
occur  within  the  zone  of  high  attenuation,  in  Figure  5  (10 
in  the  Basin  and  Range  Province  and  6  in  the  Gulf  of  Cali¬ 
fornia)  ,  and  (ii)  those  that  occur  outside  the  zone  (8  in 
the  California  Border  Province,  2  in  western  Baja,  California 
and  1  near  Denver) .  The  earthquakes  in  group  (ii)  cluster 
very  tightly  about  a  straight  line.  With  one  exception, 
earthquakes  in  group  (i)  lie  well  above  that  line.  From  the 
fitting  of  straight  lines,  using  a  least-squared-error 
criterion,  to  the  two  populations,  it  may  be  concluded  that 
in  the  mfa  range  4.5  to  5.5  earthquakes  of  a  given  surface 
wave  magnitude  in  the  Basin  and  Range-Gulf  of  California 
region  show  apparent  body  wave  magnitudes  0.3  to  0.4  less 
than  do  comparable  earthquakes  in  adjacent  areas  (principally 
the  San  Andreas  region  of  California) .  This  conclusion  receives 
support  from  the  more  recent  measurements  of  Savino  et  al. 

(1971) ;  see  Figure  1  of  their  paper.  That  this  difference 

in  the  M_  -  m,  curve  is  due  to  lateral  variations  of  Q  1 
S  o 

in  the  upper  mantle  is  strongly  suggested,  though  not  proven. 
Differences  in  average  depth  of  focus  or  source  volume  for 
earthquakes  of  the  two  areas  is  another  possible  explanation. 

The  typically  shallow  depth  of  earthquakes  in  both  regions, 
however,  and  the  finding  by  Wyss  &  Brune  (1968)  that  source 
dimensions  of  earthquakes  in  the  Nevada-Arizona  region  are 
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similar  to  or  smaller  than  those  on  the  San  Andreas  fault 
system  are  compelling  arguments  that  the  latter  explanation 
is  not  valid. 

Thus  the  attenuative  properties  of  a  medium  can  signi¬ 
ficantly  alter  the  M  -  m.  relationship,  the  criterion 
most  used  to  discriminate  underground  nuclear  explosions 
from  earthquakes.  The  anomalous  Mg  -  m^  pattern  in  western 
North  America  is  at  least  in  part  due  to  greater-than-average 
attenuation  of  P  waves  in  the  upper  mantle  of  that  region 
(cf.  Ward  &  Toksfiz  1971).  (The  single  Basin  and  Range  earth¬ 
quake  of  Mg  =  4.4  shown  in  Figure  6  with  an  unusually  large 
m^  serves  to  emphasize,  though,  that  attenuation  is  not  the 
only  parameter  controlling  the  Mg  -  m^  relation  for  a 
region.)  Even  more  importantly,  by  considering  in  Figure  6 
only  those  earthquakes  and  explosions  which  occur  in  the  same 
tectonic  region  (NTS  is  in  the  Basin  and  Range  province)  a 
more  effective  separation  of  explosion  and  earthquake  popu¬ 
lations  is  achieved  than  if  all  earthquakes  from  western 
North  America  are  lumped  into  a  single  category. 

Knowledge  of  the  frequency  dependence  of  Q  in  the  earth, 
particularly  within  the  low-Q  asthenosphere,  is  exceedingly 
important  if  body-wave  spectra  are  to  be  correctly  inter¬ 
preted  in  terms  of  source  parameters.  Seismological  evidence 
currently  supports  the  hypothesis  that  the  asthenosphere, 
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especially  in  tectonic  and  oceanic  areas,  is  partially 
melted  (Anderson  &  Sammis  1970)  .  Walsh's  (1969)  model  for 
attenuation  in  partially  melted  rock  has  been  successfully 
applied  both  to  laboratory  experiments  and  to  seismic-wave 
absorption  in  the  asthenosphere ;  a  thorough  discussion  is 
given  in  Solomon  (1972) .  The  important  predictions  of 
Walsh's  model  are  that  losses  in  compression  are  negligible 
compared  with  losses  in  shear,  verified  for  the  earth  by 
the  work  of  Anderson  et  al.  (1965)  ,  Solomon  &  Toksflz  (  1970) , 
and  others;  and  that  in  shear  behaves  as  for  a  relaxa¬ 

tion  process,  illustrated  in  Figure  7.  Both  elastic  modulus 
and  Q  are  frequency  dependent;  there  is  a  peak  in  the  atten¬ 
uation  at  an  angular  frequency  w  equal  to  1/t,  where  x 
is  a  characteristic  time  determined  in  Walsh's  model  by  such 
properties  as  the  viscosity  and  volume  concentration  of  the 
melt  and  geometry  of  the  melt-solid  interfaces.  At  frequencies 
much  less  than  1/x ,  Q  ^  is  proportional  to  frequency;  at 
frequencies  much  greater  than  1/x ,  Q  ^  is  inversely  propor¬ 
tional  to  frequency.  Note  that  the  attenuation  peak  in  a 
partial  melt  can  be  no  sharper  than  that  shown  in  Figure  7, 
according  to  Walsh's  model.  In  a  more  realistic  parameteri¬ 
zation  of  partially  melted  rock,  the  peak  would  be  "smeared 
out"  due  to  a  superposition  of  several  such  relaxations. 

An  idealized  shear-wave  spectrum  from  an  earthquake 
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source  (Aki  1967;  Brune  1970)  is  shown  in  Figure  8,  together 

with  the  spectra  expected  after  propagation  at  vertical 

incidence  through  two  simple  models  of  an  attenuating  upper 

mantle.  In  one,  Qy  =  100  independent  of  frequency;  in  the 

second,  Qy  is  given  by  Figure  7,  with  2ttt  =  10  sec  and 

0  .  =20.  In  both  models,  the  attenuating  zone  is  300  km 

mm 

thick.  At  a  frequency  of  about  1  Hz,  is  identical  in 

the  two  models. 

The  two  attenuated  spectra  are  quite  different  both  in 
absolute  amplitude  and  in  spectral  shape.  If  these  two 
spectra  were  not  corrected  for  absorption,  then  the  measured 
"corner  frequencies, "  diagnostic  of  source  dimensions  (Aki 
1967;  Brune  1970),  might  differ  by  a  factor  of  3  or  more. 

The  fall-off  of  spectral  amplitude  with  frequency  at  fre¬ 
quencies  greater  than  the  corner  value  is  related  to  frac¬ 
tional  stress  drop  in  Brune 's  (1970)  model.  Figure  8  demon¬ 
strates  the  impossibility  of  extracting  such  information  if 
the  attenuation  is  not  well  known.  We  should  note  that  the 
amplitude  changes  produced  by  the  above  two  simple  models 
for  upper-mantle  attenuation  are  something  of  a  maximum:  if 
P  waves  rather  than  S  waves  were  employed,  if  near-field 
rather  than  teleseismic  observations  were  used  to  calculate 
spectra,  or  if  the  recording  site  were  on  a  shield  instead 
of  a  tectonic  region,  then  the  effects  of  attenuation  on 
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spectral  amplitudes  would  be  less  severe. 

A  partially  melted  asthenosphere  has  an  attenuation 
spectrum  more  complicated  than  the  single  relaxation  peak 
of  Figure  7 .  In  Figure  9  is  shown  a  model  of  ,  a  function 
of  both  frequency  and  depth,  in  the  upper  mantle  of  western 
United  States  (Solomon  1972)  .  There  are  at  least  three 
distinct  relaxation  peaks,  each  operative  over  a  somewhat 
different  depth  interval.  In  Figure  10  is  given  the  ampli¬ 
tude  spectrum  of  s.  shear  wave  that  has  propagated  vertically 
through  the  upper  mantle  of  western  U.S.  Also  shown  are 
the  assumed  source  spectrum  and  the  attenuated  spectrum  for 
a  constant-Q  asthenosphere,  both  taken  from  Figure  8.  As 
noted  before,  the  corner  frequencies  and  the  overall  shape 
of  observed  body-wave  spectra  are  quite  strongly  dependent 
on  the  frequency  dependence  of  Q  in  the  upper  mantle. 

Workers  interested  in  source  spectra  generally  apply 
some  correction  for  attenuation  to  the  measured  amplitudes. 

It  is  of  interest  to  inquire  what  effect  applying  the  wrong 
Q  correction  would  have  on  the  interpretation  of  spectral 
shape.  In  Figure  10  is  shown  the  source  spectrum  that  would 
have  been  calculated  if  the  spectrum  measured  in  western 
U.S.  had  been  "corrected"  by  using  the  constant-Q  model. 

Such  a  spectrum  "corrected"  to  the  source  has  a  spurious 
relative  maximum  and  a  precipitous  decrease  in  amplitude 
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wit  h  increasing  frequency  at  frequencies  in  excess  of  the 
"peak”  value.  If  the  amplitudes  of  Figure  10  had  been 
measured  only  at  frequencies  greater  than  .1  Hz  or  so, 
the  "peak"  would  appear  to  be  an  absolute  maximum.  Inasmuch 
as  Wyss  et  al.  (1971)  have  proposed  that  such  a  peak  in  body- 
wave  amplitude  spectrum  is  diagnostic  of  an  explosion  source, 
the  need  for  properly  understanding  the  nature  of  Q  in  the 
earth  before  correcting  spectra  for  attenuation  is  clearly 
demonstrated.  It  should  be  emphasized,  however,  that  the 
above  discussion  in  no  way  alters  the  conclusions  of  Wyss 
et  al.  (1971).  Their  findings  were  based  on  the  recordings 
at  a  single  fixed  station  of  waves  from  explosions  and  earth¬ 
quakes  located  within  a  single,  small  region.  Further,  the 
spectral  peaks  reported  by  Wyss  et  al .  were  present  in  the 
data  before  any  attenuation  correction.  Note  that  a  relaxa¬ 
tion  peak  in  Q-1  can  never  produce  a  relative  maximum  in  an 
amplitude  spectrum  that  is,  prior  to  any  attenuation,  flat 
or  a  monotonically  decreasing  function  of  frequency.  This  is 

because  Q-1  never  decreases  with  increasing  frequency  f 

-1 


faster  than  as  f 


in  the  relaxation  model. 
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5 .  Conclusions 


Seismic-wave  attenuation  in  the  earth  act.s  as  a  com¬ 
plicated  low-pass  filter  for  signals  from  all  seismic  sources. 
The  important  considerations  discussed  above  that  pertain 
to  the  problems  of  describing  the  seismic  source  in  general 
and  of  discriminating  earthquakes  from  underground  nuclear 
explosions  in  particular  may  be  summarized  as  follows: 

1.  Q  in  the  earth  is  a  function  of  space  and  frequency. 

2.  The  effect  of  attenuation  on  amplitudes  is  not 
pronounced  if  all  or  nearly  all  of  the  energy  of  a  wave 
propagates  in  the  lithosphere  below  a  few  km  depth. 
Examples  include  P  and  S  waves  recorded  within  about 
1000  km  of  a  shallow  source  and  Love  and  Rayleigh 
waves  in  the  per.iod  range  30  to  about  15  sec. 

3.  If  a  major  fraction  of  the  wave  energy  propagates 
in  the  low-Q  asthenosphere,  then  proper  interpretation 
of  source  parameters  requires  consideration  of  the 
lateral  variation  and  frequency  dependence  of  Q. 


4.  Lateral  variation  of  Q  affects  only  mildly  the 
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estimation  of  surface-wave  magnitude.  After  propa¬ 
gating  3000  km  through  a  tectonic  region,  as  con¬ 
trasted  with  a  stable  continental  region,  the  surface- 
wave  magnitude  would  be  negligibly  different  at  20  sec 
and  .1  to  .3  magnitude  units  less  at  periods  greater 
than  30  sec. 


5.  The  shape  of  surface-wave  amplitude  spectra  is 
governed  to  some  extent  by  attenuation.  The  effect 
is  most  likely  to  increase  the  estimate  of  focal 
depth  from  spectral  shape  of  waves  from  explosions  if 
multi-azimuth  observations  are  not  available. 

6.  Body-wave  magnitude  is  noticeably  affected  by  the 
lateral  variation  of  Q.  The  magnitude  can  be  lower 
by  at  least  several  tenths  of  a  magnitude  unit  in 
tectonically  active  as  opposed  to  stable  regions. 

An  implication  is  that  any  comparison  of  body-wave 
amplitudes  from  earthquakes  and  explosions  should  be 
restricted  to  events  from  the  same  tectonic  province. 

7.  Body-wave  spectral  shape  is  strongly  influenced  by 
the  frequency  dependence  of  Q.  In  particular,  estima¬ 
tion  op  source  dimensions  or  fractional  stress  drop  may 
be  in  error  or  a  spurious  spectral  peak  may  be  introduced 
by  improperly  accounting  for  attenuation. 
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FIGURE  CAPTIONS 


Fractional  potential  energy  per  km  depth  stored 
per  cycle  for  Rayleigh  waves  of  various  periods 
T.  The  computational  technique  used  to  determine 
these  curves  is  from  Harkrider  (1964) .  The  velocity- 
density  model  is  the  shield  model  of  Harkrider  and 
Anderson  (1966) .  The  approximate  boundary  between 
high-Q  lithosphere  and  low-Q  asthenosphere  (shaded) 
is  also  illustrated. 

Love-wave  attenuation,  Q„ -1,  versus  frequency,  f, 

Ju  ' 

(semi-log  scale) .  The  data  of  Anderson  et  al.  (1965) 
and  Tsai  &  Aki  (1969)  are  for  primarily  oceanic  paths. 
Values  for  western  U.S.  are  from  Solomon  (1972) . 

The  surface  wave  path  for  east-central  U.S.  is  be¬ 
tween  Rapid  City,  South  Dakota,  and  Atlanta,  Georgia. 
Error  bars,  showing  standard  deviations  for  selected 
determinations,  indicate  the  precision  of  the  meas¬ 
urements  . 

Rayleigh-wave  attenuation,  QR  ^ ,  versus  frequency#  f, 
(semi-log  scale) .  Wave  paths,  error  bars,  and  sources 
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of  data  are  as  in  Figure  2. 

Fig.  4.  Rayleigh-wave  vertical  displacement  amplitude 

spectrum  from  a  nuclear  explosion  buried  .5  km  in 
tuff  (Tsai  &  Aki  1971)  and  the  spectra  expected 
after  propagation  of  such  a  wave  through  western 
and  east-central  United  States  (W.U.S.  and  E.U.S. , 
respectively) .  A  path  length  of  3000  km  is  assumed 
and  geometrical  spreading  is  ignored. 


Fig.  5.  Lateral  variation  of  the  differential  attenuation 

6t*  of  long-period  P  waves  in  the  United  States 
(from  Solomon  &  Toks&z,  1970,  with  slight  revisions). 
6t*,  in  seconds,  is  the  average  of  attenuation 
measurements  for  two  deep  earthquakes  in  South  Amer¬ 
ica.  Dashed  lines  approximate  the  6t*  =  0  contours. 


Fig.  6. 


Surface-wave  magnitude,  Mg,  versus  body-wave  mag¬ 
nitude,  m,  ,  for  events  in  southwestern  North  America, 
b 

Mg  and  m^  were  measured  at  stations  in  the  Canadian 
Network  by  Basham  (1969) .  Earthquakes  represented 


by  solid  symbols  are  located  in  regions  of  negative 
6t*  in  Figure  5;  those  given  by  open  symbols  are 
located  in  regions  of  high  attenuation  (positive  6 t* ) 
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Fig.  7. 


Fig.  8 


Fig .  9 


in  Figure  6.  Straight  lines  are  (lowest  line) 

Basham's  fit  to  nuclear  explosion  data  and  (upper 
two  lines)  least-squared  error  fits  to  the  two 
earthquake  populations. 

Elastic  modulus  M  and  attenuation  Q  1  for  a 
relaxation  process.  My  is  the  unrelaxed,  or  high- 
frequency,  modulus;  Qmin  is  the  minimum  value  of 
Q;  u  is  the  angular  frequency;  and  t  is  a  char¬ 
acteristic  time.  A  relaxation  occurs  at  wt  =  1. 

Shear-wave  spectral  amplitude  from  an  idealized 
earthquake  source  (Aki  1967;  Brune  1970)  and  the 
spectra  expected  after  propagation  of  such  a  wave 
at  vertical  incidence  through  two  simple  models  of 
upper  mantle  Q:  (i)  Q  =  100  in  a  layer  300  km  thick, 
and  is  infinite  elsewhere;  (ii)  Q  is  given  by  the 
curve  of  Figure  7,  with  2ttt  =  10  sec  and  Qmj_n  -  ' 

in  a  layer  300  km  thick,  and  is  infinite  elsewhere. 

Q  ^  in  western  United  States  (Solomon  1972) .  Qy 
is  independent  of  frequency  f  above  60  km,  and 
depends  on  frequency  as  shown  in  the  inserts  (semi¬ 
log  scale)  for  the  two  layers  60  to  160  km  deep  and 
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and  160  to  350  km  deep.  100/C^  is  plotted  as  a 
function  of  depth  z  for  two  specific  frequencies. 

Fig.  10.  Shear-wave  spectral  amplitude  from  an  idealized 
earthquake  source  (Aki  1967;  Brune  1970)  and  the 
spectra  expected  after  propagation  of  such  a  wave 
at  vertical  incidence  through  the  Q  model  of  Figure 
9  ("Western  U.S.")  and  the  constant-Q  model  from 
Figure  8  ("Q  =  100").  The  curve  labeled  '"corrected' 
to  source"  is  the  source  spectrum  that  would  have 
been  obtained  if  the  western  U.S.  curve  had  been 
corrected  for  attenuation  by  using  the  constant-Q 
model . 
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4.3  Stress  Distribution  Beneath  Island  Arcs  by  Albert 

T.  Smith,  Jr.  and  M.  Nafi  Toksttz 

Summary 

A  numerical  simulation  of  the  descending  lithosphere 
under  island  arcs  is  developed  in  order  to  describe  stress 
distribution  derived  from  earthquakes.  The  lithospheric 
slab  and  the  adjacent  mantle  are  simulated  as  an  elasto- 
static  problem  with  moduli  corresponding  to  their  effective 
rheologies.  Assuming  similar  composition  for  the  lithosphere 
and  mantle  and  given  the  temperature  regime  for  the  mantle, 
gravitational  body  forces  generated  by  thermal  volume  con¬ 
traction  of  the  colder  slab  and  by  elevation  of  phase  boundaries 
are  included.  In  addition,  simulated  convection  is  applied  to 
the  slab  in  order  to  distinguish  its  role  as  a  driving  mechanism. 
The  regional  applicability  of  gravitational  sinking  is  then 
clarified  using  a  simple  numerical  integration  of  the 
equilibrium  equation  and  comparing  this  to  available  information 
for  select  island  arcs.  These  analogues  suggest  a  viable 
hypothesis  for  the  dynamics  of  the  subjecting  lithosphere. 

First,  gravitational  sinking  induced  by  thermal  density 
anomalies  can  explain  the  directions  of  the  principle  stress 
within  the  descending  slab  at  island  arcs  without  recourse 
to  convective  effects  in  the  adjacent  mantle.  The  rheologies 
of  the  mantle  and  slab  represent  instead  the  dominant  effect. 
Second,  under  these  assumptions  the  convergence  rate  of  the 
island-arc  system  constitutes  a  major  factor  influencing  the 
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stress  distribution  by  means  of  the  gravitational  body  forces 
and  the  mantle's  resistance  or  support.  The  models  indicate 
that  a  fine  balance  exists  between  these  forces.  Furthermore , 
if  large  deviatoric  stresses  are  necessary  for  intermediate 
and  deep  earthquakes,  the  elevation  of  phase  transformation 
boundaries  are  a  necessary  feature  of  the  simulation  for  a 
mantle  and  slab  with  similar  composition.  The  post-spinel 
phase  change  at  the  650  km  depth  suggests  itself  as  a  possible 
factor  in  this  mechanism.  Finally,  a  low-strength  channel 
corresponding  to  the  seismic  low-velocity  zone  and  a  stronger 
but  constant-strength  mantle  yields  reasonable  agreement  be¬ 
tween  the  simulation  and  the  regional  stress  patterns  derived 
from  earthquakes. 
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INTRODUCTION 

New  global  tectonics  implies  both  spreading  of  the  sea 
floor  and  consumption  of  crustal  plates  into  the  mantle.  Vari¬ 
ous  geological  and  geophysical  observations  confirm  the  motion 
of  these  plates.  Furthermore,  they  impose  some  constraints 
upon  the  general  motion,  both  past  and  present.  Yet  the  mechan¬ 
ism  producing  these  motions  is  a  topic  of  considerable  debate. 
Gravitational  sinking  of  the  lithosphere,  elevation  of  the  mid¬ 
ocean  ridges,  and  viscous  drag  have  each  been  suggested  as 
the  mechanism  coupling  thermal  convection  to  the  lithosphere. 

One  means  of  limiting  the  possible  mechanisms  causing  the  motion 
is  to  compute  models  of  the  stresses  within  a  downgoing  slab, 
incorporating  various  driving  mechanisms  and  relevant  body 
forces.  Comparing  the  stress  patterns,  seismicity,  energy, 
and  focal  mechanisms  of  earthquakes  yield  constraints  upon 
the  type  of  mechanism  producing  the  plate  motion.  This  tech¬ 
nique  has  been  used  for  just  this  purpose  in  this  paper. 

Studies  of  the  source  mechanisms  of  shallow,  intermediate  and 
deep  focus  earthquakes  have  been  conducted  in  most  of  the  sub- 
duction  regions.  A  comprehensive  summary  of  the  results  are 
given  by  Isacks  &  Molnar  (1971) .  Theoretical  calculations  of 
the  stress  distributions  associated  with  a  descending  lithosphere 
have  not  kept  up  with  the  observational  data.  Using  an  analy¬ 
tical  solution  for  the  slab's  thermal  regime  (McKenzie  1969) 
approached  the  problem  of  the  stresses  within  the  slab  and 
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surrounding  mantle.  Griggs  (1971)  related  the  seismic  effects 
to  gravitational  forces  due  to  density  anomalies  arising  from 
temperature  differences  and  phase  changes  in  an  infinitely 
long  slab. 

In  this  paper  we  utilize  a  numerical  scheme  to  approach 
the  stress  problem  in  the  descending  lithosphere.  We  model 
the  descending  slab  and  the  surrounding  mantle  as  an  elasto- 
static  problem,  incorporating  body  forces  due  to  thermal  den¬ 
sity  differences  between  the  colder  slab  and  surrounding  mantle. 
Symmetry  of  the  ideal  slab  geometrically  allows  a  plane  strain 
formulation;  consequently,  the  simulation  reduces  to  an  ellip¬ 
tic  boundary  value  problem  allowing  numerical  solution.  Knowing 
the  temperature  regime  within  the  slab  (Toksflz  et.  al.  1971), 
the  induced  stresses  are  then  computed  for  elastic  moduli  se¬ 
lected  to  simulate  the  viscous  and  elastic  support  provided  by 
the  surrounding  mantle.  These  models  provide  insight  and  re¬ 
strictions  upon  the  probable  driving  mechanism,  and  in  conjunc¬ 
tion  with  other  data  can  improve  our  understanding  of  the  mo¬ 
tions  of  lithospheric  plates. 

In  the  following  sections  the  theoretical  assumptions 
and  practical  considerations  required  for  these  stress  models 
are  developed  and  the  computational  technique  is  described. 

The  dependence  of  stresses  upon  the  slab  and  mantle's  rheology 
and  their  relationship  to  source  mechanisms  of  earthquakes 
are  discussed  in  the  last  section. 
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SIMULATION  OF  LITHOSPHERIC  SLABS 

The  rheology  of  the  earth  represents  a  special  problem 
for  modeling  the  downgoing  slab.  The  lithosphere  and  the 
mantle  each  exhibit  different  behavior  to  long-term  stresses: 
the  former  acts  nearly  elastically  (Walcott  1970) ,  while  the 
asthenosphere  and  probably  the  upper  mantle  approximate  a  fluid 
with  a  low  yielding  stress  (Crittenden  1967).  Consequently , 
a  model  of  the  descending  lithosphere  must  incorporate  the 
effects  of  these  regions  and  the  dynamic  interaction  of  the 
system  composed  of  the  descending  slab  and  asthenosphere.  We 
model  the  lithosphere  as  an  elastic  model  to  reduce  this  to  a 
tractable  problem.  We  review  the  assumptions  required  for 
this  formulation  in  the  following  sections. 

Properties  of  Lithosphere  and  Descending  Slab 

Numerous  data  indicate  that  the  lithosphere  behaves 
elastically.  It  possesses  a  definite  flexure  rigidity  (Walcott 
1970)  and  displays  phenomena  such  as  bending  at  island  arcs 
(Hanks  1970)  which  are  attributed  to  its  elasticity.  Utsu 
(1967)  and  Davies  &  McKenzie  (1969)  have  observed  strong 
seismic-velocity  anomalies  at  island  arcs  corresponding  to  a 
high  velocity  channel  in  the  mantle.  This  appears  to  be  a 
direct  consequence  of  the  colder  lithosphere  descending  into 
the  mantle  (Minear  &  Toksflz  1970;  Jacobs  1970).  Furthermore, 
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the  descending  lithosphere  attenuates  seismic  waves  much  less 
than  the  surrounding  mantle  (Utsu  1967;  Molnar  &  Oliver 
1969;  Kanamori  1971,  Mitronovas  &  Isacks  1971),  again  an 
indication  of  its  elasticity.  These  and  the  observation  of 
earthquakes  support  the  assertion  of  an  elastic  lithosphere 
and  its  simulation  as  an  elastic  solid. 

Yet  it  is  not  justified  to  assume  constant  elastic  moduli 
throughout  the  slab,  for  it  is  not  a  constant  temperature. 

The  thermal  gradient  separating  the  core  of  the  slab  and  the 
mantle  must  correspond  to  a  change  in  the  elastic  moduli. 
Discussion  of  the  temperature  dependence  will  be  postponed  until 
a  later  section  in  order  to  first  assess  the  behavior  of  the 
mantle. 

Simulation  of  Mein  tie  Properties 

To  treat  the  effects  of  the  surrounding  mantle  with  an 

elastic  model,  the  mantle  must  be  modeled  as  a  'soft1  elastic 

solid  capable  of  partially  supporting  the  downgoing  slab. 

The  low  yielding  strength  and  the  viscosity  of  the  astheno- 

sphere  and  the  upper  mantle  have  been  deduced  primarily  from 

uplift  data  after  glacial  loading  (Crittenden  1967;  McConnell 

1968;  Cathles  1970).  A  low-viscosity  channel  of  1020  -  1  poise 

23 

is  generally  agreed  upon,  while  viscosities  of  10  poise  are 
indicated  for  depths  of  600  to  1200  km.  Both  experimental  and 
theoretical  results  support  the  low  creep  strength  of  the  upper 
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mantle  (Carter  &  Ave ' Lallemant  1970;  Weertman  1970;  Goetz 
1971).  Seismic  attenuation  (Anderson  et  al.  1967;  Kanamori 
1968),  electrical  conductivity  (Everett  &  Hyndman  1968),  and 
geothermal  and  melting  information  (Clark  t.  Ringwood  1964) 
substantiate  the  marked  differences  bewteen  the  lithosphere 
and  the  asthenosphere  which  has  been  attributed  to  partial 
melting  (Anderson  &  Sammis  1970).  The  upper  mantle  must  then 
be  approximated  as  &  low-strength  zone  whose  effect  upon  the 
lithosphere  arises  from  both  dynamic  stresses  and  static  sup¬ 
port.  The  simulation  must  incorporate  these  two  effects  con¬ 
sistently  with  the  geophysical  observations. 

For  our  calculations ,  the  elastic  constants  are  selected 
appropriate  to  the  degree  of  support  by  the  surrounding  mantle. 

A  rough  guide  can  be  had  from  the  effective  viscosities  and 
elastic-viscoelastic  analogy  (Christensen  1971).  Suppose  one 
observes  at  an  instant  in  time  the  stress  distribution  within 
the  mantle  and  assumes  that  the  stresses  all  originate  dynamical¬ 
ly  from  viscous  drag.  As  the  slab  is  moving  at  a  constant 
velocity,  a  very  rough  assumption  of  constant  strain  rate, 
can  be  used  with  the  effective  viscosity  H  to  define  the 
resulting  stress: 


a  =  ire 


(1) 
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This  represents  the  supporting  stress  imposed  upon  the  descen¬ 
ding  lithosphere  within  the  limitations  of  a  linear  theory. 

The  elastic  moduli  of  the  mantle  region  is  adjusted  to  fit 
these  approximate  stresses  to  yield  the  corresponding  support 
of  the  downgoing  slab. 

Figs.  4  to  10  display  Young's  modulus  for  the  com¬ 

puted  models  as  a  function  of  depth.  These  models  are  selected 
*-°  illustrate  a  range  of  mantle  conditions:  from  constant 
elastic  moduli  throughout  the  mantle  and  slab  to  a  very  low- 
strength  mantle  relative  to  the  slab.  Intermediate  cases 
having  a  'soft'  asthenosphere  and  'hard*  mantle  offer  closer 
approximations  to  the  viscous  models  deduced  by  McConnell  (1968). 
The  relative  importance  of  these  factors  can  be  judged  from 
the  resulting  models. 

Temperature  Dependence  of  Elastic  Moduli 

One  crucial  factor  remains  for  realistic  simulation  of 
the  slab:  the  elastic  constants  are  temperature  dependent. 

The  thermal  effects  upon  elastic  moduli  for  seismic  frequencies 
are  found  to  be  -  -  1.5  x  10“4/#C  for  olivine  (Chung 

1971),  where  E  is  an  effective  Young's  modulus.  But  for  long¬ 
term  stresses  this  dependence  is  not  realistic)  a  more  appro¬ 
priate  relationship  uses  the  rate  of  diffusion  creep  at  high 
temperatures.  For  this  form  of  creep  an  exponential  dependence 
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of  strain  rate  upon  temperature  is  implied  by  theoretical 
and  experimental  evidence  (Weertman  1970;  Carter  6 
Ave'Lallemant  1970) : 

e  =  f (o)exp(-Q/RT) 


where  Q  is  the  activation  energy  for  diffusion,  R  is  the 
gas  constant,  T  is  the  absolute  temperature,  and  f(o)  is  a 
function  that  depends  primarily  on  the  stress  but  that  may 
contain  other  variables.  For  constant  load  time  At  and 
stress  o  ,  it  follows  that 

E  ■  a/u  3  At°fftTexp(0/RT)  ‘  Eoexp<Q/RT)  (3> 

or  to  first  order  the  dependence  is 

“  .  - -fly  •  AT  (4> 

E  RT^ 
o 

Using  approximate  values  for  the  parameters,  a  very  rough 
upper  bound  results  for  the  descending  slab  (Weertman  1970) : 


I*  si  “  2,2  x  10  2  °c  1 


(5) 
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when  Q  =  100  kcal.#  T  =  1500  °C,  and  E  is  the  value  of 
Young's  modulus.  Figs.  5  to  7  give  models  incorporating 
a  temperature  dependence  for  the  descending  slab  lying 
within  the  bounds  of  static  (creep)  and  the  high-frequency 
limits.  These  models  provide  a  means  of  illustrating  the 
marked  effect  that  the  slab's  material  strength  has  upon 
the  resulting  stress  distribution;  one  expects  the  real 
earth  to  lie  within  this  range. 


COMPUTATION  OF  STRESSES 


Simplifying  the  rheologies  of  the  slab  and  mantle  to 
effective  elastic  moduli  allows  solving  an  elastostatic 
problem  instead  of  a  difficult  viscoelastic  problem.  With 
this  in  mind  the  representation  of  the  descending  litho¬ 
sphere  reduces  to  a  plane  strain,  elastostatic  problem 
adopting  moduli  corresponding  to  the  mantle's  support  and 
using  body  forces  based  upon  the  slab's  density  anomalies. 
Assuming  similar  composition,  the  body  forces  originate 
from  buoyancy  due  to  the  temperature  difference  between  the 
slab  and  the  adjacent  mantle.  A  region  within  the  colder 
slab  is  denser  relative  to  the  surrounding  mantle;  thus,  it 
has  a  tendency  to  sink,  thereby  inducing  stresses  within  the 
slab  and  mantle.  These  lateral  density  variations  are 
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computed  using  the  temperature  regime  obtained  by  ToksOz 
et  al.  (1971)  for  a  slab  with  a  spreading  rate  of  8  cm/yr. 

The  value  of  volume  expansion,  o  ,  is  derived  from  an 
average  of  Birch  &  Verhoogen's  (ToksOz  et  al.  1971) . 

a  =  exp (3. 58  -  0.0072z)  (6 

where  z  is  depth  in  km.  and  ct  has  units  of  10  C  . 

For  average  density  pQ  of  the  slab,  the  vertical  body 
force  Fy  due  to  negative  buoyancy  becomes 

Py  *  gp0“4T  ( 

where  g  is  the  acceleration  of  gravity,  and  AT  is  the 
horizontal  temperature  contrast  between  the  mantle  and  the 
slab.  Two  separate  values  of  the  heat  conduction  have  been 
used  for  the  temperature  calculations  illustrated  in  Fig.  1 
and  2.  One  set  of  models  uses  MacDonald's  formulation  of  the 
radiative  term  (ToksOz  et  al.  1971),  while  the  second  set 
has  a  reduced  radiative  term  derived  from  Schatz's  (1971) 
experimental  investigations  for  olivine.  In  the  latter  case 
the  radiative  conductivity  is  given  by 

K_  =  0  for  T  <  500°K,  and 
R 

K_  -  5.5xl0~6 (T-500.)  for  T>500°K 
R 


(8) 
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In  addition  the  effect  of  the  phase  change  is  included  in 
the  body  force  as  a  separate  problem.  The  depth  and  density 
changes  of  the  olivine-spinel  transition  are  estimated 
using  Ringwood  &  Major's  data  (1970),  and  the  elevation 
of  the  boundary  is  derived  from  the  thermal  regime  in 
accordance  with  Toksflz ' s  calculations .  With  the  inclusion 
of  these  parameters,  the  problem  reduces  to  an  elliptic 
boundary  value  problem  with  variable  coefficients. 

A  finite-difference  scheme  with  a  20  km  grid  spacing 
has  been  used  to  solve  this  problem  for  an  840  by  680  km 
region  enclosing  the  slab  and  adjacent  mantle.  The  digital 
computer  program  has  been  developed  for  on-line  solution 
of  generalized  elliptic  boundary  value  problems  by  C. 

Tillman  (1971)  using  a  nine-point,  integral  formulation  of 
the  difference  equations  and  successive  point  over-relaxation 
to  obtain  a  solution.  Utilization  of  this  program  has 
allowed  solution  of  the  equilibrium  and  constitutive  equation 
for  an  inhomogeneous,  isotropic  body  with  infinitesimal 
deformation  (Nowinski  &  Turski  1953) : 

aa8,8  =  “Fa(xl'x2}  (a'6  =  1,2) 

vE(xx,x2)  E(x1,x2) 

°a8  =  (1  +  v)  (1  -  2v)  6a86  +  2(1  +  v)  (uct,B  +  Ua,8} 


(9) 
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u  s  a  component  of  displacement 

U 

v  =  Poisson  ratio  E(x^,X2)  =  Young's  modulus 

Reducing  these  to  the  divergence  form  in  terms  of  displacements, 
one  has  a  set.  of  coupled  equations: 


3xl  (1  +  v)  (1  -  2\>) 


3  r  E(x,y) 
3yl2(l  +  v) 


u  )] 
y  ,x 


(10) 


3  .  E  (x ,  y) 
3xl2(l  +  v) 


+ 


9_r  E(x,y) 

3y LT1  +  v) (1  -  2v)U x,x 


+ 


(1  -  v)E(x,y) 
(1  +  v)(l-2v) 


“y.y1 


“Fy(x,y) 


where  Fa  =  a  component  of  body  force, 


3u . 

_ l 

3x  . 
D 


or 


u. 


x,y 


etc. 


The  spatial  variations  of  body  forces  and  elastic  parameters 
are  introduced  and  a  solution  follows  from  the  difference 
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scheme  with  the  appropriate  boundary  conditions. 

If  the  boundaries  surrounding  the  slice  of  mantle  are 
far  enough  from  the  slab,  their  effect  upon  the  descending 
lithosphere  will  be  minimal.  This  corresponds  to  St. 
Vincent's  principle  in  elastostatics.  The  two  extreme  cases 
are  a  rigid  boundary  (no  displacements)  and  a  free  surface 
(no  tractions) .  For  the  side  boundaries,  the  computations 
indicate  minimal  effect  of  the  boundaries  upon  the  stress 
distribution  within  the  descending  slab.  Our  previous 
arguments  for  a  relatively  soft  mantle  suggest  small 
tractions  at  the  side  boundaries;  hence,  these  have  been 
selected  for  our  computations.  The  top  face,  corresponding 
to  the  earth's  surface,  is  also  specified  as  free,  whereas 
the  bottom  face  is  made  rigid.  This  lower  boundary  is 
necessary  to  satisfy  the  requirements  for  a  static  problem: 
the  sum  of  all  couples  and  forces  must  be  zero  for  the  region. 
Again  the  boundary  does  not  seriously  effect  the  stresses 
within  the  slab.  In  addition,  some  models  have  the  edges 
of  the  lithosphere  rigid  to  simulate  the  junction  of  the 
descending  slab  and  the  crustal  plate.  Yet  due  to  the 
limitations  of  these  boundary  conditions,  the  lack  of 
information  for  the  lithosphere  as  it  bends  into  the  mantle, 
and  the  characteristics  of  the  plate  boundaries,  these 
models  cannot  hope  to  indicate  the  stresses  for  depths  less 
than  about  100  km,  and  that  is  not  their  intent.  Rather,  the 
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intermediate  and  deep  regions  of  the  downgoing  slab  are 
within  the  realistic  domain  of  the  models. 

Stability  and  convergence  of  the  solution  is  ensured  by 
both  test  problems  and  on-line  inspection  of  convergence. 

The  program  has  successfully  solved  a  variety  of  problems 
including  potential  flow,  bending  of  plates,  and  others 
(Tillman  1971) .  This  series  has  confirmed  the  utility 
and  stability  of  the  difference  scheme.  In  addition, 
inhomogeneous,  plane  strain  elastic  solutions  confirm  the 
stability  for  problems  such  as  the  descending  slab  (appendix) . 
The  20  km  grid  spacing  used  for  the  slab  is  adequate  for  the 
broad  features  encountered  in  the  thermal  regime  and  for  the 
desired  resolution  of  the  stress  distribution.  Finally, 
on-line  monitoring  and  selection  of  the  relaxation  parameters 
provides  for  optimum  convergence  and  continual  inspection 
of  stability  (Mitchell  1969) . 


RESULTS  OF  SIMULATION 

Theoretical  Results 

A  set  of  eight  models  simulating  different  rheologies 
and  conditions  are  summarized  in  Table  1  together  with  the 
important  parameters  and  results.  Figs.  3  through  10 
display  these  models  including  the  mantle  and  slab's  elastic 
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moduli.  each  figure  depicting  the  variation  of  a  particular 
parameter.  The  differences  between  the  two  thermal  regimes 
in  Fig  1  and  2  are  not  the  controlling  factor  in  these 
models;  the  reduced  radiative  conductivity  increases  the 
gravitational  body  force  and  density  gradient,  yet  this  is 
a  minute  effect  compared  to  the  changes  in  mantle  support 
and  slab  rheology.  The  mantle  and  slab  properties  control 
the  behavior  of  the  stress  pattern:  the  maximum  shear  stress 
concentrates  within  the  colder,  harder  core  for  a  slab  with 
temperature  dependent  elastic  parameters.  Decreasing  the 
mantle's  total  and  distributed  support  alters  the  principle 
stress  from  compression  downgoing  to  tension  along  the  trend 
of  the  slab.  These  relationships  will  be  clarified  in  the 
following  discussion  of  the  ndividual  models. 

Comparison  of  Models 

A  simulation  is  first  necessary  for  a  control  case, 
specifically,  a  model  having  identical  elastic  moduli  within 
the  mantle  and  slab.  Model  1  in  Fig.  3  represents  such  a 
case.  For  this  case  the  maximum  shear  stress  occurs  at  the 
lower  edge  of  the  slab  and  within  the  adjacent  mantle,  whose 
stresses  have  not  been  contoured.  The  orientation  and  the 
maximum  stress  reflect  the  distribution  of  body  forces;  at 
the  lower  edge  and  within  the  adjacent  mantle,  the  com- 
pressional  axis  is  vertical  as  one  expects  from  simple  static 
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consideration.  The  distribution  of  body  forces  also  intro¬ 
duces  downdip  tension  for  shallow  depths.  Given  this  example, 
the  succeeding  models  illustrate  the  marked  effects 
produced  by  different  effective  rheologies  within  the  slab 

and  mantle. 

Models  2  and  3  in  Fig.  4  represent  slabs  with  constant 
properties  and  different  surrounding  mantle  conditions.  The 
slab  can  be  thought  of  as  a  bending  plate  with  different 
degrees  of  support  along  its  length.  In  model  2  the  mantle 
differs  little  from  the  slab,  and  only  a  small  reduction  in 
strength  occurs  in  the  asthenosphere.  The  slab  is 
relatively  immobile  in  the  deep  mantle,  whereas  the  astheno¬ 
sphere  allows  'sagging*  analogous  to  a  plate  held  at  its 
edges:  compression  along  the  top  face  and  tension  along  the 

bottom.  If  the  asthenosphere  is  'softer'  as  in  model  3, 
the  slab  almost  hangs  as  though  suspended  by  the  lithosphere 
aud  only  partially  supported  at  its  base.  Compression  then 
occurs  at  the  bottom  face  of  the  slab  and  tension  along  its 
upper  surface.  Deceasing  the  mantle's  support  also  increases 
the  stress  guided  through  the  slab.  Yet  the  shear  stresses 
are  concentrated  at  the  edges  of  the  slab  as  one  expects  for 
a  bending  plate  but  not  necessarily  for  the  descending 
lithosphere. 

An  alternate  set  of  models  is  desired  which  concentrates 
the  shear  stress  towards  the  slab's  center;  models  4  and  5  in 
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Fig.  5  show  such  a  transition.  The  slab's  strength  is 
strongly  dependent  upon  temperature  according  to  equation  (3 
with  a  pressure  term  included: 

Q  +  pV 

E  *  exp( — — — )  ,  or  for  Vac  =  0, 

Eslab  =  5,0  X  lol°exP(-0-36xlo~2AT) 

+  Emantle 

where  p  is  the  pressure  in  kilobars,  and  Vgc  is  the 
activation  volume  for  diffusion  creep  (Weertman  1970) .  The 
stresses  are  now  strongly  focused  within  the  colder  regions 
of  the  slab  as  expected  from  diffusion  creep.  In  addition, 
two  extreme  cases  of  mantle  support  are  illustrated:  one 
having  the  tensional  or  minimal  principle  axis  following 
the  trend  of  the  slab  along  its  full  length,  the  other  with 
the  compressional  axis  paralleling  the  slab.  The  downgoing 
slab  in  model  5  hangs  from  the  lithosphere.  With  increasing 
strength  and  support  by  the  mantle,  the  slab  can  be  placed 
under  compression  as  in  model  4.  A  balance  then  exists 
between  the  mantle's  support,  the  slab's  strength,  and  the 
distribution  of  density  anomalies  within  the  slab.  The 
next  set  of  models  clarifies  this  balance  between  the  mantle 
and  the  slab's  rheology. 
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Sup  pose  the  strength  of  the  slab  is  only  weakly  dependent 
upon  temperature  as  given  by: 


Eslab 


E 


mantle 


+  5 . OxlO10  [1 


(100.  +  0.895p)AT] 
2(Tmantle)  xl°  (12) 


Model  6  in  Fig.  6  simulates  such  a  condition  for  a  mantle 
with  a  low-strength  channel  similar  to  model  4.  Its  stress 
distribution  is  an  intermediate  case  between  the  downdip 
compression  of  model  4  and  the  bending  stresses  that  are 
characteristic  of  model  2;  however,  the  temperature  dependence 
of  the  elastic  moduli  is  insufficient  to  focus  the  stress 
towards  the  center  of  the  slab.  As  in  model  4  the  narrow 
asthenosphere  prevents  tension  along  the  slab's  full  length. 
Model  7  is  similar  to  model  6  but  has  no  temperature 
dependence  for  the  descending  lithosphere.  As  expected,  its 
maximum  shear  stress  concentrates  at  the  edges  of  the  slab; 
however,  the  soft  asthenosphere  introduces  tensional 
features  at  intermediate  depths.  The  models  suggest  that  a 
critical  value  for  the  temperature  dependence  may  be 
necessary  in  order  to  focus  the  stress  within  the  slab. 

On  the  other  hand,  increasing  the  contrast  of  the  elastic 
moduli  within  the  slab  strongly  focuses  the  stress  and  can 
induce  tension  at  intermediate  depths.  Fig.  7  illustrates 
this  case  with  model  8  for  a  temperature  dependence  given 
by 
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E  ~  E  -i  exp[  -(1.2x10  +  5.5x10  p)  AT]  (13) 

slab  mantle 

A  very  soft  asthenosphere  combined  with  a  hard  central  core 
for  the  slab  produce  high  shear  stresses  with  the  tensional 
axis  following  the  slab  at  intermediate  depths  and  with 
downdip  compression  for  greater  depths. 

Three  parameters  are  resolved  that  best  described  these 
categories.  The  first  measures  the  ratio  of  mantle  strength 
at  two  depths,  140  and  440  km,  thereby  resolving  the 
differential  support  of  the  asthenosphere  and  the  deeper 
mantle.  Yet  this  is  incomplete  without  a  measure  of  the 
contrast  between  the  mantle  and  the  slab,  as  given  by  the 
ratio  of  their  elastic  constants  at  460  km  depth.  The 
final  parameter  measures  the  temperature  dependence  of  the 
slab's  elastic  moduli.  It  is  the  ratio  of  moduli  at  its 
edge  compared  to  its  center  at  160  km  depth.  Fig.  8  gives 
a  three  dimensional  perspective  of  these  models  in  terms 
of  the  parameters.  The  categories  are  representative  of 
the  model  space  spanned  by  this  series  of  computations. 

Thus,  regions  are  blocked  off  depending  upon  the  models 
locations  and  characteristics.  The  boundaries  are  merely 
to  illustrate  the  regional  dependence  and  should  only  be 
interpreted  as  guidelines. 
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Phase  Changes  and  Convection 

Fig.  9  considers  the  elevation  effect  of  the  olivine- 
spinel  phase  change  within  the  slab  (Schubert  &  Turcotte 
1971) .  Using  Ringwood  &  Major's  (1970)  estimation  of 
the  slope  for  the  phase  transformation  boundary,  30  bar/  C. 
the  density  anomalies  within  the  slab  are  computed  using 
Fig.  1  when  the  transition  occurs  over  50  km  (ToksOz  et  al. 

1971)  .  The  problem  then  lends  itself  to  superposition: 
calculating  the  solution  independently  and  summing  the 
stresses  for  the  separate  models.  The  transition  induces 
local  stress  concentrations  directly  above  and  below  the 
perturbed  depth.  The  slab  tends  to  be  under  tension  at  or 
above  350  km  and  in  compression  below  this  depth.  Consequently, 
if  the  slab  is  originally  under  compression  at  350  km,  the 
transition  reduces  the  maximum  shear  stress  above  that 
region  while  increasing  it  below  the  phase  change. 

In  the  same  manner  the  convection  problem  in  Fig.  10  is 
treated  as  a  slab  with  boundary  tractions  along  its  faces.  The 
right-hand  problem  does  not  include  the  mantle's  support; 
instead,  the  slab  acts  as  a  free  plate  with  applied  tractions. 
The  figure  illustrates  an  example  using  50  bars  resistance 
along  the  slab's  upper  surface  and  no  tractions  along  its 
bottom  face.  This  roughly  simulates  convection  along  the 
bottom  face  and  viscous  drag  on  the  top.  As  expected  for  a 
plate,  bending  effects  are  superimposed  upon  the  dominate 
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horizontal  principle  stress.  The  left-hand  modal  with 
mantle  support  incorporates  a  shear  couple  along  the  upper, 
mantle-slab  boundary  to  simulate  drag.  The  stress  orientations 
are  now  more  complicated  due  to  the  mantle  support.  Yet  the 
models  reveal  that  such  convection  generally  increases  com¬ 
pression  along  the  trend  of  the  slab. 

Comparison  to  Island  Arcs 

Isacks  &  Molnar  (1971)  have  summarized  the  observed 
stress  patterns  for  island-arc  regions  and  their  relationship 
to  the  regional  tectonics  and  gravitational  sinking  of  the 
descending  lithosphere.  Fig.  11  gives  their  interpretation 
of  the  principle  stresses  derived  from  focal-mechanism 
solutions.  Three  crucial  regions  are  apparent  in  the  stress 
pattern  and  the  numerical  models:  at  shallow  depths  the 
tensional  principal  stress,  T,  parallels  the  slab's  orientation; 
intermediate  orientations,  B,  or  down-dip  extension  occurs 
at  intermediate  depths;  and  down -dip  compression,  P, 
dominates  for  greater  depths.  The  following  sections  will 
show  the  correspondence  between  the  computed  models  and  the 
regional  stress  distribution,  and  in  addition,  will  corroborate 
the  sinking  lithosphere  hypothesis.  Regions  have  been 
selected  as  representative  of  island-arc  systems  and  for 
adequate  knowledge  of  the  arc,  a  prerequisite  for  discriminating 
the  crucial  elements  of  the  mechanism.  Unless  referenced  to 
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the  contrary,  the  seismicity  and  the  focal-mechanism  solutions 
are  contained  in  Isacks  4  Molnar  (1971)  . 

Honshu 

The  dipping  slab  beneath  North  Honshu  and  the  Sea  of 
Japan  closely  resembles  model  4:  a  30  degree  dipping  slab 
under  compression  for  depths  greater  than  80  km.  The 
seismic  activity  is  continuous  and  does  not  extend  much  beyond 
600  km  in  depth.  Concentrating  upon  the  relatively  un¬ 
contorted  central  region,  compression  is  preferred  for  the 
majority  of  solutions.  Closer  inspection  of  the  data  reveals 
a  slight  rotation  towards  down-dip  tension  of  the  P  axis 
with  decreasing  depth,  yet  nothing  that  can  be  termed 
definitive.  Kanamori's  (1971)  solution  of  the  Sannku  earth¬ 
quake  of  1933  has  down-dip  extension  and  is  located  near  the 
trench  axis.  Only  for  shallow  depths,  then,  is  the  slab 
under  tension;  the  bulk  remains  under  down-dip  compression 
as  attested  by  the  focal-mechanism  solutions. 

Comparing  these  orientations  to  those  of  model  4  in 
Fig.  5,  the  principle  stress  is  also  down-dip  compression  for 
intermediate  and  greater  depths.  The  simulation  implies  that 
a  narrow,  low-strength  asthenosphere  can  exist  if  sufficient 
support  is  provided  by  the  metasphere.  In  addition,  the 
model  suggests  rotation  of  the  stress  orientations  for  inter¬ 
mediate  depths.  The  nature  of  these  rotations  cannot  be 
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determined  on  the  basis  of  a  two-dimensional  simulation  of 
the  real  earth.  Moreover,  Kanamori's  (1971)  tensional 
solution  for  the  Sanriku  earthquake  demonstrates  that  the 
simulation  is  not  fully  satisfactory;  the  theoretical  model 
is  still  compressional  for  shallow  depths.  Thio  effect 
has  its  origin  from  either  of  two  related  factors:  the 
slab's  core  is  insufficiently  'hard',  or  its  corrolary,  the 
mantle  is  supporting  an  excessive  fraction  of  the  slab’s 
body  force.  The  major  factor,  however,  could  be  the  rigid 
boundary  at  the  edge  of  the  lithosphere.  It,  in  effect, 
simulates  the  lithosphere  pushing  the  descending  slab. 
Consequently,  a  simulation  having  parameters  intermediate 
to  models  4  and  8  is  necessary  for  extension  at  shallow 

depths . 

Tonga 

The  marked  similarity  between  Tonga  and  Honshu  in  Fig.  11 
dictates  a  corresponding  interpretation  of  the  models;  both 
involve  support  by  the  mesosphere  and  a  soft  asthenosphere 
as  suggested  by  Isacks  &  Molnar  (1969).  Both  convergence 
rates  are  approximately  9  cm  per  year  (Le  Pichon  1968) 
implying  greater  density  anomalies  and  depth  before  equilibrium 
(Minear  &  Toksttz  1970) ,  higher  strength  within  the  colder 
slab,  alterations  of  phase  boundaries,  and  perhaps  greater 
resistance  or  drag  by  the  mantle.  Their  interrelationship 
clarifies  itself  when  additional  regions  are  examined. 
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Kermadec 

The  Kermadec  arc,  lying  just  south  of  Tonga,  represents 
a  marked  change  in  the  stress  distribution.  The  inclined 
zone  parallels  a  slab  dipping  60  degrees  west  and  extends 
to  500  or  550  km  depth  (Sykes  1966) .  Limited  focal  plane 
solutions  show  down-dip  extension  at  230  km  depth  and 
compression  for  350  km.  The  region  corresponds,  then,  to 
model  8  in  Fig.  7,  if  the  mantle  support  is  further  reduced 
to  allow  extension  at  230  km  depth.  Thus  in  Fig.  8  the 
regions  are  evolving  towards  decreasing  mantle  support. 
Moreover,  the  convergence  rate  at  the  Kermadec  arc  is 
markedly  less  than  the  preceeding  two  cases,  5.4  cm/yr  as 
opposed  to  9  cm/yr . 

South  America 

Peru  and  Northern  Chilean  regions  carry  this  deduction 
one  step  further:  a  gap  in  seismicity  exists  between 
intermediate  depth,  extensional  earthquake  mechanisms  and 
compressional  mechanisms  at  600  km.  Briefly,  the  shallow  ox 
intermediate  zone  is  characterized  by  down-dip  tension  for  a 
shallow  dipping  slab;  however,  no  conclusive  seismic  activity 
occurs  between  200  or  300  km  and  500  km  for  Peru  and  Chile, 
respectively.  Between  500  and  600  km  the  seismic  activity 
dramatically  increases,  and  the  mechanisms  indicate  down-dip 
compression.  Whether  the  slab  is  continuous  through  the 
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seismicity  gap  remains  a  conjecture.  In  Chile  high-frequency 
shear  waves  pass  through  the  gap,  while  the  results  are 
inconclusive  for  Peru  and  western  Brazil  (Molnar  &  Oliver 
1969;  Sacks  1969).  Either  the  zone  of  attenuation  in  the 
mantle  occurs  above  the  gap,  or  the  descending  lithosphere 
is  continuous  to  depths  of  650  km. 

Assuming  the  slab  is  continuous,  the  seismicity  and 
focal  mechanisms  could  follow  from  the  dynamics  of  gravi¬ 
tational  sinking.  The  North  Chilean  slab  extends  200  km 
deeper  than  Kermadec,  but  the  convergence  rate  is  similar, 

5.2  cm/yr  (Le  Pichon  1968)  .  The  additional  length  of  slab 
has  two  opposing  effects;  the  total  body  force  is  increased 
by  the  added  segment  and  by  the  contribution  of  any  deep 
phase  changes  within  the  slab.  Yet  deeper  penetration  of  the 
slab  and  variations  in  the  mantle's  resistance  or  support 
opposes  the  increased  body  force.  An  extrapolation  of  model 
8,  which  has  tension  at  intermediate  depths,  could  produce 
such  a  gap  given  the  proper  combination  of  mantle  and  slab 
support,  body  forces,  and  perhaps  elevation  of  phase-transformation 
boundaries.  The  additional  length  of  slab  increases  the 
total  body  force  and  offsets  any  mantle  support  gained  by 
greater  penetration.  A  balance  must  exist  between  the  body 
force,  mantle  resistance,  and  depth  of  penetration  for  the 
occurrence  of  this  pattern.  For  Kermadec  and  Chilean  slabs, 
the  resistance  or  support  at  each  depth  is  held  constant  by 
similar  convergence  rates;  consequently,  the  effects  caused 
by  the  additional  length  of  slab  are  present  given  that 
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their  preceding  history  is  inconsequential.  The  latter 
assumption  is  questionable.  Otherwise,  if  the  slab  is 
continuous  and  if  the  convergence  rate  has  been  reasonably 
stable  over  the  past  10  million  years,  it  is  an  adequate 
proposition. 

Other  evidence  from  South  America  suggests  that  the 
slab  behaves  as  a  stress  guide  for  both  shallow  and  deep 
earthquakes.  Carr  et  al.  (1971)  reported  evidence  of 
lithospheric  bending  for  shallow  depths:  tensional  earth¬ 
quakes  near  the  surface  and  down-dip  compression  along  the 
bottom  face.  This  cannot  be  attributed  to  thrusting  of 
one  lithospheric  plate  under  another.  The  Kurile  arc  also 
shows  compression  for  shallow  depths  at  the  lower  face.  For 
greater  depths  Fukao  (1971)  has  obtained  the  solution  for  a 
deep- focus  earthquake  (577  km)  in  western  Brazil  including 
the  fault  plane  parameters.  It  verifies  shear  faulting  for 
deep  earthquakes  and  indicates  that  the  P  axis  is  nearly 
vertical,  corresponding  to  the  trend  of  the  slab.  The  stress 
drop,  300  to  460  bars,  implies  that  the  compressional 
principle  stress  is  greater  than  600  to  920  bars.  Moreover, 
the  faulted  area,  510  to  580  km  ,  the  average  dislocation, 

265  to  355  cm,  and  the  location  of  subsequent  earthquakes 
suggests  large-scale  faulting  over  nearly  the  entire  thickness 
of  the  descending  lithosphere.  These  observations  are 
consistent  with  gravitational  sinking  as  the  mechanism  inducing 
earthquakes  within  the  descending  lithosphere. 


-184- 


Middle  America  and  the  Aleutians 

The  Middle  American  arc  corresponds  to  a  short  slab 
dipping  about  40  to  60  degrees  at  depths  greater  than  85  km. 
Focal-mechanism  solutions  are  consistent  with  down-dip 
extension  along  its  200  to  300  km  length.  Thus,  the  arc 
resembles  model  5  with  minimal  support  from  the  surrounding 
mantle. 

The  Aleutian  arc  also  conforms  to  a  shallow-slab 
representation  with  a  dip  of  approximately  60  degrees  (Davies 
&  McKenzie  1969)  and  with  seismicity  generally  less  than 
200  to  300  km  deep.  Unlike  the  Middle  American  arc, 
compressional  solutions  are  evident  for  intermediate  depths, 
while  down-dip  tension  occurs  for  shallow  depths  adjacent  to 
the  trench  (Stauder  1968a, b) .  The  arc  could  illustrate  a 
truncated  version  of  model  8,  similar  perhaps  to  model  7; 
however,  the  factors  that  determine  the  individuality  of  the 
Aleutians  and  Middle  American  arcs  remain  obscure.  Conceivably, 
the  distinctive  stress  distributions  result  from  the 
increased  convergence  rate  of  the  Aleutian  arc,  5.2  cm/yr 
(Le  Pichon  1968)  compared  to  roughly  3.5  cm/yr  for  Middle 
America  (Molnar  &  Sykes  1969) ,  thereby  augmenting  the  thermal 
body  forces  and  the  mantle  resistance  and  placing  the  slab 
under  partial  compression.  These  convergence  rates  are, 
however,  questionable. 
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Kurile -Kamchatka 

The  Kurile-Kamchatka  arc,  dipping  at  45  degrees,  poses 
a  crucial  validation  of  these  mechanisms.  Contrary  to  the 
other  island  arcs,  the  seismicity  suggests  the  maximum 
depth  of  the  zone  is  probably  less  than  600  km  and  perhaps 
nearer  500  km.  Yet  the  convergence  rate  is  approximately 
7.5  cm/yr  (Le  Pichon  1969).  Based  upon  the  proposed 
mechanism,  one  would  expect  a  distribution  similar  to 
Kermadec  with  its  500  km  slab  but  with  addition?!  body  forces 
and  mantle  resistance  generated  by  the  faster  convergence 
rate.  Indeed,  the  effect  is  apparent  in  the  stress  dis¬ 
tribution:  a  minimum  in  activity,  and  not  an  aseismic  gap, 
is  reported  between  about  200  and  400  km.  The  shallow 
events  above  the  gap  exhibit  down-dip  extension,  whereas 
for  events  located  beneath  the  gap  the  stress  orientation 
is  down-dip  compression.  Moreover,  rotation  of  the  principle 
axis  occurs  for  events  near  the  periphery  of  each  region 
just  as  in  model  8.  These  characteristics  place  the  arc  as 
an  intermediate  case  between  Kermadec  and  Tonga *  and 
substantiate  the  idea  that  the  descent  velocity  of  the  litho¬ 
sphere  plays  a  major  role  in  the  dynamics  and  manifests  itself 
in  the  seismicity  and  stress  orientations  by  means  of  the 
body  force  and  mantle  resistance. 


-186- 


DISCUSSION  AND  IMPLICATIONS  -  A  MODEL  OF  DESCENDING  LITHOSPHERE 

The  results  discussed  in  the  previous  section  suggest  a 
simple,  semi-empirical  analogue  to  explain  the  regional 
stress  patterns  observed  for  island  arcs.  When  the  earthquake 
characteristics,  the  limiting  depth  of  down-dip  tension  and 
compression,  are  plotted  in  terms  of  the  convergence  rate,  a 
relationship  appears  between  the  seismicity,  focal  mechanisms, 
and  the  convergence  rate.  Fig.  12  illustrates  this  pattern 
in  addition  to  a  theoretical  model.  Beginning  with  New 
Zealand  and  slow  convergence  rates,  a  large  seismicity  gap 
occurs  between  down-dip  tensional  solutions  and  deep  com- 
pressional  solutions.  For  faster  convergence  rates  with  a 
deep  slab,  the  zone  narrows  and  then  disappears  for  the  Kurile 
island  arc,  yet  a  minimum  in  seismicity  persists  between 
200  and  400  km  in  addition  to  fluctuations  in  the  stress 
orientations.  Tonga  also  displays  a  minimum  within  this 
range.  The  second  figure  indicates  an  analogous  pattern  for 
steeply-dipping  slabs.  Indeed,  the  relationship  can  be 
deduced  from  gravitational  sinking  of  the  descending  litho¬ 
sphere  if  a  continuous  slab  is  presupposed. 

A  simulation  of  gravitational  sinking  suggests  the 
dependence  of  the  stress  distribution  upon  the  body  force  and 
the  mantle  resistance  via  the  convergence  rate.  Postulating 
the  dependence  of  body  forces  and  mantle  support  upon  the 
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convergence  rate  of  the  island  arc,  and  assuming  the  depth 
dependence  for  the  support  and  phase  changes,  a  numerical 
integration  of  the  forces  along  the  slab  yields  the  stress 
within  the  slab  as  a  function  of  convergence  rate,  length 
of  slab,  and  dip.  The  stresses  are  only  those  along  the 
direction  of  the  slab  and  indicate  whether  tension  or  com¬ 
pression  prevail.  Fig.  14  portrays  a  numerical  solution 
given  the  following  assumptions:  first,  the  body  force  due 
to  thermal  contraction  is  constant  along  the  slab  and  varies 
as  the  square  root  of  the  rate.  It  has  been  normalized  to 
yield  140  dynes/cm3  at  8  cm/yr;  the  precise  value  does  not 
affect  the  solution.  The  rate  function  roughly  corresponds 
to  the  models  of  Minear  &  Toksfiz  (1970)  and  equals  the 
rate  of  conduction  for  an  infinite  body  wich  initial 
temperature  distribution  (Carslaw  &  Jaeger  1958) .  Alterations 
of  the  phase  boundaries  are  also  given  this  rate  dependence 
when  normalized  to  8  cm/yr  boundaries  calculated  by  Toksttz 
et  al.  (1971) .  Viscous  drag  is  now  assumed  for  the  support 
or  resistance;  thus,  resistance  is  linearly  dependent  upon 
the  rate.  In  addition,  the  resistance  at  its  end  is  2  to  5 
times  greater  than  at  the  faces.  Increasing  tension  with 
decreasing  convergence  rate  as  in  Fig.  12  suggests  that 
resistance  must  decrease  faster  than  the  total  body  force  if 
the  slab  is  continuous  through  the  aseismic  gap.  The 
dependence,  controlled  by  mantle  rheology,  -is  then  selected 
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to  yield  the  desired  distribution.  Finally,  the  resistance 
is  normalized  at  10  cm/yr  to  obtain  compression  throughout 
the  slab.  To  compute  steeply  dipping  slabs,  only  the 
component  of  body  force  along  the  dip  is  altered,  while 
other  factors  such  as  greater  temperature  contrasts  remain 
constant.  Restrictions  upon  the  length  of  the  downgoing 
slab  and  the  convergence  rate  have  not  been  imposed  and  are 
not  necessary  in  this  analogue  (Luyendyk  1970) .  Table  2 
summarizes  the  parameters  and  solution  for  one  successful 

model . 

Certain  observations  are  in  order  regarding  the 
characteristics  of  successful  models:  if  stresses  greater 
than  1  kilobar  are  deemed  necessary  (Wyss  1970;  Fukao  1971) 
and  if  compositional  differences  are  insignificant,  the 
additional  body  force  provided  by  elevated  phase  boundaries 
is  a  prerequisite.  Alone,  body  forces  from  thermal 
contraction  are  insufficient  to  eauso  the  seismicity  and 
the  large  apparent  stresses  observed,  particularly  at  great 
depths  in  South  America.  Perhaps  the  deep  seismicity, 
maximum  at  600  km  reflects  a  change  in  the  material  propertie 
of  the  slab  as  could  occur  in  the  post-spinel  phase  change. 
Only  the  depth  of  the  phase  change  must  decrease  within  the 
slab.  A  positive  slope  of  dP/dT  for  the  post-spinel  phase 
transformation  would  imply  both  higher  body  forces  and  a 
change  in  physical  properties  at  600  km. 
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The  analogues  also  employ  a  low-strength  asthenosphere 
extending  to  220  km;  however,  little  or  no  increase  in 
resistance  is  demanded  for  greater  <epths.  This  is  in  accord 
with  a  recent  solution  of  glacial  uplift  data  (Cathles 
1971) ,  but  directly  contradicts  models  of  gravitational 
sinking  relying  upon  monotonic  increasing  strength  with 
depth  along  the  length  of  the  slab  (Isacks  &  Molnar  1969, 
1971)  . 

Assuming  a  continuous  slab  for  all  seismic  zones 
entails  no  presuppositions  for  the  mechanism  of  separation. 
Instead,  the  aseismic  gap  stems  from  its  low  stress. 

Neither  hypothesis  subjects  the  deeper  portions  to  stress; 
both  must  rely  upon  forces  adjacent  to  and  within  the  deep 
slab  as  the  mechanism  inducing  deep  earthquakes  in  South 

America  and  similar  locales. 

Fig.  13  further  corroborates  the  model's  relationship 
to  seismicity.  A  distinct  correlation  appears  between  the 
regional  seismicity  and  the  stress  levels  within  the 
hypothetical  slab.  Only  the  relative  shapes  of  the  curves 
are  significant;  the  actual  magnitude  of  stress  within  the 
slab  depends  upon  the  precise  model.  Underthrusting  of 
lithospheric  plates  also  contributes  to  shallow  seismicity, 
and  variations  in  the  material  properties  of  the  slab  may 
be  important  for  all  depths. 

Finally,  the  models  readily  explain  the  pattern  for 
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steeply  dipping  slabs:  the  greater  down-dip  component  of 
the  body  force  pulls  the  zone  of  tension  further  down  the 
slab.  Yet  constant  resistance  along  the  length  of  the 
slab  and  localized  body  forces  from  phase  transitions 
allow  compression  even  for  nhort  slabs.  This  emphasizes 
the  importance  of  these  phase  transformations  upon  the 
proposed  slab  dynamics. 


CONCLUSIONS 

Summary 

Gravitational  sinking  induced  by  thermal  density  anomalies 
can  explain  the  directions  of  the  principle  stress  within  the 
descending  slab  at  island  arcs  without  recourse  to  signif¬ 
icant  convective  effects  in  the  adjacent  mantle.  Rather, 
the  rheologies  of  the  mantle  and  slab  are  the  dominant 
factors  influencing  the  stress  pattern.  Earthquake  source 
mechanisms  together  with  the  models  indicate  that  the 
temperature  dependence  of  the  slab's  material  properties  is 
significant;  bending  effects  that  are  characteristic  of  a 
plate  are  only  observed  for  shallow  depths.  Instead,  the 
stress  is  transmitted  predominately  down-dip  within  the  slab. 

A  low-strength  asthenosphere  or  channel  below  the  lithosphere 
is  also  consistent  with  the  models  and  with  the  regional 
earthquake  mechanisms.  These  are  in  agreement  with  Isack3  & 


Molnar  (1971) . 

Numerical  integration  for  a  simple  model  containing 
the  relevant  forces  and  their  assumptions  suggests  a  hypo¬ 
thesis  for  the  dynamics  of  subduction.  An  essential  element 
is  the  dependence  upon  the  velocity  of  the  converging  litho¬ 
spheric  plates;  it  represents  a  major  influence  upon  the 
dynamics  of  descent.  The  thermal  body  forces,  boundaries  of 
phase  changes,  and  elastic  constants  are  manifestations  of 
the  convergence  rate.  Yet  a  crucial  dependence  results 
between  the  mantle's  support  or  resistance  and  the  velocity 
in  order  to  account  for  the  regional  distribution  of  earth¬ 
quakes:  the  resistance  increases  with  faster  convergence 

rates.  The  viscous  end-resistance  may  also  be  more  important 
than  hitherto  acknowledged.  Finally,  increasing  strength 
with  depth  in  the  mantle  is  not  a  prerequisite  for 
simulating  the  regional  stress  distribution.  A  low-strength 
channel  corresponding  to  the  seismic  low-velocity  zone  and 
a  stronger  but  constant  strength  mantle  yields  reasonable 
agreement  between  the  models  and  the  regional  stress  patterns 
derived  from  earthquakes.  If  the  mantle  and  slab  have 
significantly  different  compositions,  the  simulations  are  no 
longer  valid. 

In  addition,  the  contribution  of  phase  changes  is 
predicted  as  a  crucial  element  of  the  hypothesis  given  the 
large  apparent  stresses  observed  for  intermediate  and  deep 
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earthquakes  and  the  initial  assumptions  for  the  models.  The 
hypothesis  is  then  consistent  with  a  form  of  shear  faulting 
as  the  dominant  source  mechanism  for  these  earthquakes. 

The  simple  models  suggest  further  that  the  post-spinel  phase 
change  at  650  km  could  provide  a  mechanism  needed  to  induce 
deep  earthquakes,  given  a  positive  slope  for  dP/dT  at  the 
phase  boundary.  Otherwise,  the  high  stresses  proposed  for 
deep  earthquakes  are  probably  inconsistent  with  the 
hypothesis . 

The  computational  results  of  this  paper  rely  upon  a 
variety  of  theoretical  and  experimental  inferences  that  require 
further  research. 
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Appendix 

FORMULATION  OF  FINITE-DIFFERENCE  EQUATIONS 

The  finite-difference  approximation  for  equation  (10) 
witnin  the  digital  program  (Tillman  1971)  uses  the  self¬ 
adjunct  form  of  the  elliptic  equations  together  with  an 
integral  formulation  for  the  difference  equations  (Varga 
1962) .  A  complete  development  may  be  found  in  the 
references ;  however ,  we  will  review  the  formulation  for 
the  elastostatic  problem. 

The  self -adjunct,  elliptic  boundary  value  problem  has 
the  form 


E<i£  +  !#  +  +  3y<e3x  +  + 


(A-l) 


+  Gu  =  H 


where  x  and  jr  are  the  independent  variables,  u  —  (u^  ,Uj *  •  •  •  Ujj) 
is  the  set  of  unknown  functions  to  be  determined,  and 
A  =  (a^)  ,  B  =  (bkl)  ,  . . . ,  G  =  (3^-l)  and  H  =  (hj^)  ,  k,  1  =  1,  1, 
. .  .N,  are  specified  matrices  whose  components  may  vary  with 
x  and  The  integral  formulation  of  equation  (A-l)  reduces 
to  an  area  integral  over  any  subregion  (Varga  1962;  Tillman 
1971)  : 
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On  application  of  Gauss's  theorem,  this  becomes 
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where  the  boundary  Z  of  subregion  R  is  traced  in  the 
clockwise  direction. 

For  the  elastostatic  equations,  the  formulation 
reduces  to  a  set  of  difference  equations  for  a  uniform 
orthogonal  lattice  with  grid  lines  parallel  to  the  x,;£ 
axes  and  separated  by  a  distance  h.  The  interior  node, 

Nq,  is  surrounded  by  8  adjacent  nodes,  each  numbered 
consecutively  starting  at  the  x  axis  and  working  clockwise. 
In  this  case  the  set  of  interior  difference  equations 
at  an  interior  node  NQ  becomes : 
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where  <5(z)  is  a  difference  approximation  to  <z.  Thus, 

<S  (g^)  is  a  vector  approximation  to  gjj-  in  terms  of 
the  approximations  to  the  dependent  variables  u  =  (u^ ,  Uj) 
at  node  and  its  neighbors  ,  N2***»  Ng*  The 

summation  is  along  the  boundary  E  of  the  subregion.  This 
subregion  is  a  square  centered  upon  the  interior  node  N0 
and  with  side  h.  The  subregion  may  be  further  subdivided 
into  8  equal  right  triangles,  tk,  all  having  a  vertex 
in  common  at  the  interior  node  Nq.  s^,...Sg  denotes  the 
boundary  of  subregion  starting  at  the  x  axis  and  continuing 
clockwise,  and  Ayk,  ^xj,  is  hhe  length  of  the  segment  sk 
in  the  x  and  y.  directions.  Sk  indicates  that  the  matrix 
or  difference  approximation  is  evaluated  at  the  midpoint 
of  segment  sk,  and  Tk  requires  evaluation  at  the  centroid 
of  the  triangular  subregion  tk.  The  evaluation  of  the 
governing  system  of  matrices  A,  C, . . . ,  H  and  the  approximation 
6 (g~)  ,  etc.,  are  accomplished  by  a  bilinear  interpolation 
from  the  values  at  the  surrounding  nodes.  Thus, 


-196- 


evaluation  for  the  midpoint  of  segment  requires  a 

linear  interpolation  using  nodes  NQ,  N^,  and  . 

The  formulation  results  then  in  a  nine-point  difference 
approximation  to  equation  (10)  for  the  displacement  field. 
The  matrices  for  equation  (A-4)  become 


\ 
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H  = 


-Fy(x,y) 


(A-5) 


E(x,y)  represents  the  spatially  varying  Young's  modulus,  and 
F  (x,y)  is  the  gravitational  body  force.  Using  the 
finite-difference  approximation  to  the  coupled  equations, 
a  solution  may  be  obtained  using  successive  over-or-under- 

relaxation. 

The  finite-difference  equations  must  be  represented 
in  matrix  form 


where  A  and  B  are  a  coefficient  matrix  and  inhomogeneity 
vector,  respectively,  produced  by  means  of  the  integration 
technique.  The  vector  of  approximations  to  the  nodal 
values  of  the  dependent  variables  u,  is  determined  by 
successive  over-relaxation.  For  a  problem  with  N  dependent 
variables,  each  node  has  a  set  of  N  equations  which  are 
ordered  successively  in  (A-6)  according  to  the  grid  numbering. 
Correspondingly,  the  unknown  vector  U  of  the  system  is 
formed  by  stacking  sets  of  nodal  approximations  Up  -  <u1'u2)p 
end-on-end  in  the  order  of  the  numbering.  Thus  for  a 
lattice  with  M  nodes,  the  vector  U  of  (A-6)  would  have  MN 
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components  while  A  would  be  a  square,  MN  by  MN,  coefficient 
matrix,  and  B  an  inhomogeneity  vector  with  MN  components. 

The  solution  is  then  obtained  using  the  method  of  successive 
point  over-  or  under-relaxation  (Varga  1962) . 

The  solution  of  the  general  plane  elastostatic  problem 
given  by  equation  (A-4)  and  (A-5)  requires  theorems  for 
stability  and  convergence  of  coupled  elliptic  equations 
with  variable  coefficients.  Since  these  are  not  available, 
we  must  depend  upon  empirical  tests  of  convergence  (Mitchell 
1969)  . 

Mild  restrictions  are,  however,  imposed  by  the 

y 

integration  technique  (Varga  1962;  Tillman  1971)  ,: 

i 

(1)  dependent  variables  u  =  (u^,..»,  un)  must;  be 
continuous  functions  of  x  and  y; 

(2)  all  coefficients  of  the  equations  must  be 

•  ’  v 

piecewise  continuous; 

(3)  there  must  be  flux  continuity  in  any  direction 
and  at  any  point  of  the  region. 

Flux  continuity  is  actually  a  property  imposed  upon  the 
solution  u  i'/  the  integration  technique  and  guarantees  that 
proper  irharface  conditions  are  formulated  along  lines  of 
material  discontinuities  without  special  attention  (Varga 
1962) . 

Convergence  towards  a  true,  discrete  solution  using 
successive  over-relaxation  depends  upon  the  properties  of  the 
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coefficient  matrix,  which  in  turn  depends  on  the  particular 
differential  system  and  finite-difference  lattice  under 
cons ider>  t ion.  Successive  relaxation  will  converge  for 
any  relaxation  factor  w  in  the  range  0<w<2.  Tf  matrix 
A  is  symmetric,  nonsingular,  and  positive  definite,  and 
further,  convergence  will  occur  for  w  =>  1  (Gauss-Seidel 
method)  if  A  is  strict'' y  or  irreducibly  diagonally  dominant 
(Varga  1962).  The  last  condition  is  satisfied;  however, 
showing  that  matrix  A  possesses  the  other  properties  is  a 
difficult  task,  particularly  for  variable  coefficients. 

One  must  then  rely  upon  empirical  tests  to  evaluate  the 
convergence . 

Tillman  (1971)  describes  several  sample  problems 
successfully  solved  by  the  program.  These  include 
conduction  problems,  simple  plane  stress,  potential  flow, 
and  linear  and  nonlinear  plate  bending,  including 
irregular  lattice  structures.  Comparison  to  exact  solutions 
available  for  some  problems  indicates  insignificant  errors; 
the  final  answers  are  within  a  few  percent  of  the 
theoretical. 

Fig.  14  depicts  the  numerical  solution  for  a  plane 
strain  problem  with  variable  elastic  moduli.  A  sguare  with 
11  grid  points  to  the  side  is  taken  for  the  region,  and 
Young's  moduli  is  made  to  vary  along  the  y-axis.  A  rigid 
boundary  is  placed  along  the  left-hand  face;  free  sides  are 
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specified;  and  25  dynes/cm2  outward  traction  is  imposed  by 
means  of  body  forces  upon  the  right-hand  face.  Thus  the 
square  is  under  tension.  The  exact  solution  along  the 
axis  of  the  square  corresponds  within  6%  to  the  numerical 
results  for  u  -  1.5  and  after  150  iterations.  The  problem 
indicates  that  on-line  monitoring  can  provide  an  accurate 

evaluation  of  the  convergence. 

For  general  solutions ,  empirical  evaluation  must 
be  relied  upon  to  determine  convergence,  for  it  does  not 
necessarily  occur  for  all  values  of  the  relaxation  parameter. 
Yet  on-line  tests  for  different  parameters  furnish  a 
reliable  appraisal  of  the  convergence;  both  convergence  to 
a  discrete  solution  and  the  optimum  rate  may  be  evaluated. 
This  technique  has  proven  the  surest  means  of  guaranteeing 
stability  and  convergence  for  the  problems. 
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Table  I-A.  Models  of  Descending  Lithosphere  using  Elastostatic  Simulation 
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Boundary  condition  on  edge  of  lithosphere 


Table  I-B.  Stress  Distribution  in  Descending  Lithosphere  Model 
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Down-dip  stress:  P,  compression;  T ,  tension;  I,  intermediate. 

Maximum  shear  stress:  low,  <150  bar;  med,  150-300  bar;  high,  >300  bar 
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Table  II.  Numerical  Integration  of  Equilibrium  Equation  for 


Descending 

Slab:  45°  dip?  8  cm/yr; 

80  km 

thick 

Down-din  Stress* 

Depth 

Body**  Mantle 

(km) 

Force  Resistance 

slab 

(dyne/cm^) 

(bar) 

depth:  700 

560 

280 

35. 

-100. 

30. 

896 

2209 

672 

71. 

-100. 

30. 

434 

1747 

210 

106. 

-100. 

30. 

-29 

1285 

-252 

141. 

-100. 

30. 

-491 

823 

-714 

177. 

-100. 

144. 

-953 

361 

-1176 

212. 

-100 

296. 

-1273 

41 

-1496 

247 

-100. 

606. 

-1405 

-91 

-1628 

283. 

-180. 

606. 

-1148 

166 

-1371 

318. 

-260. 

606. 

-1291 

23 

353. 

-260. 

606. 

-1833 

-520 

389. 

-260. 

606. 

-2376 

-1063 

424. 

-100. 

606. 

-2919 

-1606 

459. 

-100. 

606. 

-2662 

-1349 

495. 

-100. 

606. 

-2405 

-1092 

530. 

-180. 

606. 

-2148 

-834 

565. 

-260. 

606. 

-2291 

-977 

601. 

-260. 

608. 

-2834 

636. 

-260. 

715. 

-3374 

671. 

-100. 

822. 

-3780 

707. 

-100. 

929. 

-3252 

End- re 

sis Lance  (bar) : 

-2591 

-1520 

-1514 

*  Positive  denotes  tension  (bar) 
**  Component  along  the  slab 
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Figures 

Fig.  1. 

Temperature  field  at  time  t  =  12.96  by  using 
MacDonald's  (1959)  value  of  the  radiative  conductivity. 
Adiabatic  compression,  phase  changes,  shear-heating,  and 
radioactive-heating  are  included  as  energy  sources  in  the 
model.  The  spreading  rate  is  8  cm/yr.  Shading  indicates 
zones  of  phase  changes.  For  the  elastostatic  problem, 
the  slab  is  truncated  when  its  temperature  reaches 
equilibrium  with  the  surrounding  mantle  at  560  km. 

(reproduced  from  Toksttz  et  al.  1971) . 

Fig.  2. 

Temperature  field  at  time  t  =  9.45  my  for  a  30  degree 
dipping  slab  with  Schatz's  (1971)  radiative  conductivity. 

All  other  parameters  are  the  same.  The  MacDonald  geotherm 
has  been  used  for  the  mantle.  Notice  the  larger  thermal 
gradients  and  temperature  contrasts  within  the  slab 
(provided  by  N.  Sleep  1971) . 

Fig.  3. 

Maximum  shear  stress  for  model  1:  contours  of  maximum 
shear  stress  in  bars  are  shown  in  the  slab  for  stresses  less 
than  500  bars.  The  arrows  indicate  the  direction  of 
deviatoric  compression  or  tension  (principle  stresses  minus 
hydrostatic  pressure)  for  select  nodes  of  the  finite-difference 
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solution.  In  this  case,  Young’s  modulus,  E,  is  constant 
throughout  the  region  for  both  the  mantle  and  slab,  and 
it  equals  1.5  1012  dynes/cm2.  The  gravitational  body 
forces  result  from  thermal  contraction  using  the  temper¬ 
ature  regime  in  Fig.  1. 

Figs.  4a  and  b. 

Solutions  for  the  maximum  shear  stress  for  two  different 

mantle  conditions  using  the  temperature  regime  of  Fig.  1. 

The  contours  on  the  right  are  again  those  of  maximum  shear 

stress  in  bars  within  the  slab.  Shear  stresses  greater  than 

500  bars  are  not  contoured  in  these  illustrations.  The 

stresses  within  the  mantle  are  not  contoured.  On  the  left, 

2 

E  represents  Young's  modulus  (dynes/cm  )  used  for  the 
mantle  surrounding  the  slab.  Unless  otherwise  specified, 
the  elastic  moduli  of  the  slab  are  constant  within  it  and 
equal  to  the  lithosphere's  value. 

Fig.  4a 

Model  2  depicts  sagging  of  the  slab  at  160  km  in  the 
asthenosphere:  compression  near  the  top  face,  tension  at 
the  lower  surface.  Only  a  slight  decrease  in  the  astheno- 
sphere's  strength  is  used  for  the  model. 

Fig.  4b 

In  contrast,  model  3  with  low  strength  asthenosphere 
and  increasing  strength  with  depth. 
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Figs.  5a  and  b. 

Solution  of  maximum  shear  stress  for  two  different 
mantle  supports  using  a  30  degree  dipping  slab  in  Fig.  2. 
Section  A-A*  depicts  the  temperature  dependence  of  Young's 
modulus  within  the  slab  for  both  models. 

Fig.  5a. 

Model  4  with  a  soft  asthenosphere  and  increasing  mantle 
support.  Compression  results  along  the  full  length  of  the 
slab. 

Fig.  5b. 

Model  5  with  little  mantle  support.  The  slab  hangs 
under  tension. 

Fig.  6a. 

Maximum  shear  stress  for  model  6  using  the  thermal 
regime  in  Fig.  2  and  a  slowly  varying  thermal  dependence  in 
the  slab  given  by  section  A-A'.  Sagging  effects  and  stress 
concentrations  again  result  along  the  edge  of  the  slab. 

Fig.  6b. 

Maximum  shear  stress  for  model  7  when  the  thermal 
regime  of  Fig.  1  is  truncated  at  380  km.  An  extended 
asthenosphere  is  used  for  a  slab  with  constant  elastic 
moduli.  Both  sagging  and  intermediate  tensional  stresses 
occur  within  the  slab. 


*  T\ 


-213- 


Fig.  7. 

Maximum  shear  stress  for  model  8  illustrating  a  slab 
with  strongly  varying  elastic  moduli  given  by  section 
A-A'.  The  thermal  regime  in  Fig.  1  is  used  for  the 
computations.  Large  down-dip  tensional  stresses  pre¬ 
dominate  at  intermediate  depths  and  gradually  rotate 
to  compression  for  greater  depths. 

Fig.  8. 

Three-dimensional  diagram  suimarizing  the  dependence 
of  the  models  upon  the  elastic  moduli.  The  three 
coordinates  represent  increasing  mantle  strength  with 
depth  relative  to  the  asthenosphere  (mantle/mantle) , 
increasing  temperature  dependence  of  the  slab's  elastic 
moduli  (slab/slab) ,  and  increasing  contrast  between  the 
slab  and  mantle's  elartic  moduli  (slab/mantle).  The 
points  refer  to  models  1  through  8  in  'rable  1.  The  regions 
blocked  off  are  meant  only  as  a  rough  guide  to  the 
characteristic  stress  patterns  observed  for  the  previous 
models. 

Fig.  9. 

The  maximum  shear  stress  for  model  9  with  olivine- 
spinel  phase  change  introduced  at  400  km  depth  according 
to  Ringwood  &  Major's  (1970)  data  and  using  the  thermal 
regime  in  Fig.  1.  A  density  anomaly  of  approximately 
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0.3  gm/cm3  results  at  350  km.  The  elastic  moduli  in  the 
mantle  and  slab  are  the  same  as  model  2. 

Fig.  10. 

Maximum  shear  stress  for  two  models  simulating 
convection.  Model  10  on  the  right  is  for  a  free  plate 
with  50  bars  traction  along  its  upper  face  simulating  mantle 
resistance.  The  arrows  again  indicate  direction  of 
principle  stress  at  various  nodes.  The  maximum  shear 
stress  is  less  than  100  bars  throughout  the  slab.  Model 
11  on  the  left  uses  a  shear-couple  at  the  mantle-slab 
interface  to  simulate  resistance  by  the  mantle.  The 
arrows  indicate  the  component  of  the  shear  couple  acting 
on  the  slab,  while  a  similar  component  but  in  the  opposite 
direction  acts  upon  the  mantle.  The  mantle's  elastic 
modulus  for  this  model  is  given  on  the  left  and  is  the 
same  as  model  2.  The  maximum  shear  stress  is  approximately 
70  bars  throughout  the  slab. 

Fig.  11. 

Global  summary  of  the  distribution  of  down-dip  stresses 
in  inclined  seismic  zones.  The  stress  axis  that  is 
approximately  parallel  to  the  dip  of  the  zone  is  represented 
by  an  unfilled  (open)  circle  for  compression  or  P  axis, 
and  a  filled  (solid)  circle  for  the  tensional  or  T  axis; 
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an  "X"  indicates  that  neither  the  P  nor  the  T  axis  is 
approximately  parallel  to  the  dip.  For  each  region  the 
line  represents  the  seismic  zone  in  a  vertical  section 
aligned  perpendicular  to  the  strike  of  the  zone.  The 
lines  show  approximately  the  dips  and  lengths  of  the  zone 
and  gaps  in  the  seismic  activity  as  a  function  of 
depth  (reproduced  from  Isacks  6  Molnar  1971) . 

Figs.  12a  and  b. 

Down-dip  stresses  in  the  inclined  seismic  zones  as 
a  function  of  convergence  rate  and  depth.  The  filled 
circles  indicate  the  maximum  depth  of  down-dip  tension; 
the  open  circles  show  the  shallowest  down-dip  compressional 
solution  (Isacks  &  Molnar  1971) •  If  two  of  the  same 
circles  occur  for  a  region,  they  represent  the  lower  and 
upper  bound  for  the  deepest  tensional  solution  or  the 
shallowest  compressional  event  inferred  from  the  boundaries 
of  the  aseismic  gap  in  Fig.  11.  The  Kurile  island  arc 
apparently  has  tensional  and  compressional  solutions 
scattered  throughout  the  region.  The  convergence  rates  are 
taken  from  Le  Pichon  (1968,  1970).  The  lines  across  the 
diagram  depict  theoretical  limits  when  the  down-dip  stress 
is  greater  than  1  kilobar.  The  solid  line  defines  the  lower 
bound  for  tension;  the  dashed  line  is  the  upper  bound  for 
compression;  and  stresses  generally  less  than  500  bars  occur 
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between  the  limits.  The  maximum  depth  of  the  slab  is 
given  on  the  right  for  each  theoretical  model. 

Fig,  12a. 

Limits  of  tension  and  compression  for  inclined  seismic 
zones  dipping  less  than  45  degrees.  Theoretical  models 
for  a  45  degree  dipping  slab  are  included  for  slabs 
penetrating  560  and  700  km  in  depth. 

Fig.  12b. 

Inclined  seismic  zones  and  theoretical  models  for 
steeply  dipping  slabs  (60  degrees  or  more) .  New  Zealand 
has  been  included  in  the  previous  diagram  since  the  curves 
merge  for  slowly  converging  slabs. 

Fig.  13a. 

Annual  number  of  earthquakes  per  25  km  depth  intervals 
as  a  function  of  depth  for  three  select  regions:  Tonga, 
Kurile-Kamchatka,  and  Alaska  (data  from  Sykes  1966) . 

Fig.  13b. 

Magnitude  of  down-dip  principle  stress  as  a  function 
of  depth  for  slab  models  using  numerical  integration  of 
the  equilibrium  equation.  The  gravitational  body  forces 
and  mantle  resistance  are  a  function  of  convergence  rate. 
The  models  correspond  to  the  previous  three  regions :  a 
700  km  depth  slab  at  9  cm/yr  convergence  rate  for  Tonga; 
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560  km  depth  at  8  cm/yr  for  the  Kurile  island  arc?  and 
300  km  depth  at  6  cm/yr  for  Alaska-Aleutian  slab.  The 
dips  are  shown  along  side  the  depth.  Effects  of  litho¬ 
spheric  bending  and  overthrusting  are  not  included  for 
shallow  depths.  The  relative  shapes  of  the  curves  show  a 
distinct  correlation. 


Fig.  14. 

Finite-difference  solution  for  a  plane-strain 
elastostatic  problem  with  variable  elastic  moduli  in  the 
y-direction.  The  boundary  conditions  for  the  square 
region  are  illustrated  at  the  bottom.  The  arrows  indicate 
the  25  dyne/cnf  body  force  imposed  at  the  boundary  nodes 
on  the  far  right.  The  20  cm  sides  consist  of  11  grid 
points.  The  solution  along  the  y-axis  is  shown  at  the 
top  together  with  the  value  of  Young’s  modulus.  The 
theoretical  solution  is  25  dynes/cm2  except  near  the  rigid 

boundary . 
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